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SUMMARY 


A method has been developed for two- and three-dimensional computations 
of viscous supersonic jet flows interacting with an external flow. The 
approach employs a reduced form of the Navier— Stokes equations which allows 
solution as an initial-boundary value problem in space, using an efficient 
noniterative forward marching algorithm. Numerical instability associated 
with forward marching algorithms for flows with embedded subsonic regions is 
avoided by approximation of the reduced form of the Navier— Stokes equations 
in the subsonic regions of the boundary layers. Supersonic and subsonic 
portions of the flow field are simultaneously calculated by a consistently 
split linearized block implicit computational algorithm. The results of 
computations for a series of test cases associated with supersonic jet flow 
is presented and compared with other calculations for axisymraetric cases. 
Demonstration calculations indicate that the computational technique has 
great promise as a tool for calculating a wide range of supersonic flow 
problems including jet flow. Finally, a User’s Manual is presented for the 
computer code used to perform the calculations. 


I. INTRODUCTION 


During the past two decades much effort has been expended in developing 
numerical procedures which can be used as alternatives to solving the full 
Navier-Stokes equations for certain classes of problems (Ref. 1). These 
procedures treat a reduced form of the steady state Navier-Stokes equations, 
often referred to as the 'parabolized Navier-Stokes equations,' as an initial 
boundary value problem that can be solved by spatial forward marching. The 
ability to obtain a solution by forward marching the governing equations from 
an initial streamwise location to some desired downstream location rather 
than perform a global solution of the governing equations, as is required for 
the solution of the full Navier-Stokes equations, results in a considerable 
savings of computational time. Although the amount of savings will depend on 
the problem considered, the efficiency of the solut ion procedures and 
numerous other variables, this savings is the primary motivation for the 
development of these marching procedures. 

To devise a set of governing equations suitable for the spatial forward 
marching of supersonic flows, three steps must be taken. First, a nominal 
primary flow direction must be identified. Second, a coordinate system must 
be constructed with one of its coordinate directions closely aligned with the 
primary flow direction. Third, all diffusion in the primary flow direction 
must be neglected. These steps when applied to the steady Navier-Stokes 
equations produce a set of governing equations which is well posed for the 
spatial forward marching of supersonic flows (e.g. Ref. 2). The introduction 
of no slip surfaces into a supersonic flow results in the formation of 
embedded subsonic regions adjacent to these surfaces. When the set of 
reduced equations, without further approximation, is forward marched with 
embedded subsonic regions the governing equations are not well posed and 
hence the solution procedure may become unstable. Even when the flow in this 
embedded subsonic region is approximated further and governed by what are 
essentially the boundary layer equations, an instability can still be 
encountered. This particular instability, which is often referred to as the 
branching phenomenon, has been the subject of much research (e.g. Refs. 3-8) 
and the technique used to surpress this instability is a convenient way to 
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differentiate between procedures for solving the reduced form of the 
Navier-Stokes equations for supersonic flow with embedded subsonic regions. 

In one of the earliest works in this area Garvine (Ref. 3) demonstrated 
(for a model problem) the existence of exponentially growing (divergent) 
terms in the spatial development of a solution of a reduced form of the 
Navier-Stokes equations when applied to the problem of an inviscid supersonic 
flow interacting with a viscous boundary layer. The author concluded that 
for this problem the reduced form of the Navier-Stokes equations was 
improperly set as an initial value problem, because the interaction dynamics 
contained upstream “elliptic** influence. In the model problem, if the 
upstream conditions are not precisely set as to cause the divergent terras to 
be multiplied by zero, the exponentially growing terms will cause the 
streamwise pressure gradient terms to grow exponentially large resulting in 
unrestrained acceleration or deceleration of the flow. In general it is not 
possible to pick the upstream conditions to negate the exponentially growing 
modes, hence several investigators have attempted to suppress the unstable 
(or branching) behavior by further modification of the reduced form of the 
Navier-Stokes equations for supersonic flows with embedded subsonic regions . 

Much of the early work on the solution of the reduced form of the 
Navier-Stokes equations is based on the work of Rudraan and Rubin (Ref. 4). 
Rudman and Rubin solved the equations for the hypersonic flow over slender 
bodies with sharp leading edges. Based on a order of magnitude analysis they 
demonstrated that for this class of problems the streamwise pressure gradient 
term was negligible when compared with the inertia and viscous terms of the 
streamwise momentum equation. Neglecting the streamwise pressure gradient 
term together with all streamwise diffusion results in a reduced form of the 
Navier-Stokes equations that is well posed for spatial forward marching even 
with embedded subsonic regions and branching was not observed in their 
calculations. Although this approach does yield a set of equations that is 
well posed for spatial forward marching, the assumption of negligible 
streamwise pressure gradient limits the class of flows which can be 
considered. In a later work Lubard and Helliwell (Ref. 5) proposed a method 
for preventing branching that involved explicit spatially lagged evaluation 
of the streamwise pressure gradient term. When marching from the i^ to 
the i+l st streamwise station all streamwise terms (See Fig. 1), except the 
streamwise pressure gradient, are evaluated by a backward difference. The 
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streamwise pressure gradient term is approximated by differencing the 
streamwise pressure gradient at prior known spatial locations, i.e. at 
the i-l st and i t ^ station (hence the terminology explicit evaluation). 

The above authors found that in addition to the frequently encountered 
problem of instability associated with exceeding some marching direction step 
size, a further instability is encountered when the step size is reduced 
below some limit. By examining the eigenvalues of a model set of equations 
(Ref. 9) they were able to develop a criterion for this minimum step size. 
(Numerical experimentation with their computer code demonstrated reasonable 
correlation with their criterion). Numerous flow fields have been 
successfully predicted using this method (mainly for cone flow) by the 
authors of Ref. 5 and others (Refs. 10-12), and in these cases evidently 
the restriction on the minimum marching step size was not a problem in 
allowing sufficiently accurate results to be obtained. However, the 
restriction on minimum marching step size is, in principle, not a desirable 
feature, since it does prevent arbitrary mesh refinement, and thereby the 
assurance that an accurate unique solution has been obtained. In at least 
one case (Ref. 9) this minimum step size restriction prevented the authors 
from successfuly obtaining a solution. In a later technique developed by 
Rakich, Vigneron, and Agarwal (Ref. 6) a variant of the technique of 
Lubard and Helliwell was used to prevent branching. In this particular 
variant the streamwise pressure gradient term is approximated by an implicit 
backward difference in the supersonic portion of the flow. However, in the 
subsonic region only that portion of the streamwise pressure gradient term 
that can be included without causing branching is evaluated implicitly. 

The results of a stability analysis similar to that of Ref. 9 also produces a 
restriction on the minimum allowable step size. When that portion of the 
subsonic pressure gradient which could not be evaluated implicitly is 
evaluated explicitly by a lagged technique similar to Lubard and Helliwell, 
Rakich et al noted that the scheme became unstable. Thus, in order to 
achieve stability this technique neglected the explicit portion of the 
streamwise pressure gradient term in the subsonic region, and implicitly took 
into account only that portion of the term that can be stabily computed. 
Schiff and Steger (Ref. 7) treat the subsonic streamwise pressure 
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gradient term by what the authors term either a first- or second-order 
streamwise extrapolation technique in the subsonic regions. The first order 
technique is equivalent to setting the streamwise pressure gradient term 
equal to zero in the subsonic region while the second order technique is 
equivalent to the explicit evaluation of the streamwise pressure term (as was 
done by Lubard and Helliwell). As with the two previously discussed 
techniques, these authors also report a restriction on the minimum marching 
step size that they may take and still retain a stable calculation. Lin and 
Rubin (Ref. 8) have developed a global relaxation procedure for solving the 
reduced form of the Navier-Stokes equations. This technique was primarily 
developed for application to cases where upstream influence is strong. To 
obtain the upstream influence with the reduced form of the Navier-Stokes 
equations requires a global iteration or relaxation procedure. The above 
authors do this by approximating the streamwise pressure gradient term by a 
forward difference. When marching the solution from the i 1 -^ to the i+l st 
station the pressure gradient term is evaluated in terms of the pressure at 
the i+l st and i+2 nd station (the i+l st station is the implicit station; 
all other streamwise derivatives are backward differenced between i+1 
and i). Initially the (unknown) pressure at the i+2 nd station is guessed; 
during subsequent global iterations the previously calculated value is used. 
Global iteration of the governing equations is continued until the solution 
converges. Lin and Rubin report that convergence is typically obtained in 
five to ten iterations for cases with small streamwise pressure gradients 
(cases run to date have been limited to flow over cones). The authors also 
report that there is no minimum marching step size requirement with their 
approach . 

The purpose of the present investigation is to develop an efficient 
numerical procedure for the solution of the two and three-dimensional reduced 
form of the Navier-Stokes equations for high Reynolds number internal flow. 
The study is limited to cases where the incoming flow is supersonic and the 
flow inside the internal flow device is, in the mean, supersonic. The 
existence of embedded subsonic regions adjacent to the surfaces of the 
internal flow devices is to be accounted for as part of the analysis. 

Because of the complexity of the physical processes occurring in internal 
flow devices and, especially in three dimensions, the large number of grid 
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points required (and hence computer time) to accurately resolve these 
processes, it was decided that a technique that solves the full Navier-Stokes 
equations would be used only if no suitable alternative could be found. 

The physics of supersonic internal flow devices is characterized by the 
formation of shock waves, the growth of boundary layers, and the interaction 
of these phenomena. Many of the above internal flow phenomena are turbulent 
and have associated with them large streamwise pressure gradients, e.g., a 
high Reynolds number incident shock wave-boundary layer interaction. It is 
to be expected that in regions of such an interaction one would desire to 
take a small streamwise marching step to accurately resolve the phenomenon. 
In particular, it might prove necessary to resolve the turbulent boundary 
layer viscous sublayer (large cross flows can occur in this region) and take 
marching steps of this order. Reviewing the existing methods for Refs. 4-8 
causes one to be concerned that techniques having such a minimum step size 
might not permit sufficient resolution of the large gradients expected. The 
method of Rudman and Rubin presumes that the streamwise pressure gradient is 
small in comparison with the inertia and viscous terms. For the cases they 
considered this is a valid assumption. This is not the case for a shock 
wave-boundary layer interaction in moderately supersonic flow. The methods 
of Refs. 5-7 all make an attempt to consider the effect of the streamwise 
pressure gradient in the embedded subsonic regions. However, they give only 
an approximate treatment to this possibly dominant term and all of those 
methods have a minimum marching step size limitation which, in many cases of 
interest in this study, may not either allow for an accurate or in some cases 
even a minimally acceptable solution. Subsequently, it will be demonstrated 
that for a case with a large streamwise pressure gradient the minimum step 
size size limitation of the order of magnitude found in Ref. 9 was 
insufficient to accurately resolve the phenomenon. Here we seek a 
noniterative approach with a consequent reduction in computer cost relative 
to either the global iteration approach to solving the reduced form of the 
Navier-Stokes equations or solution of the full Navier-Stokes equation. 
Further as a prerequisite, we require that there exist no numerical 
limitation on the minimum marching step and it is desired to keep to a 
minimum any approximation to the streamwise pressure gradient term. 
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In view of the above, it was decided to develop a numerical procedure 
for the solution of the reduced form of the Navier-Stokes equations, with 
special emphasis to be placed upon application to internal flow devices. The 
remainder of this report will describe that effort. It will consist of 
(1) a discussion of the analysis used in the study, (2) a discussion of the 
solution of the governing equations, (3) the results of a series of test 
cases run to demonstrate the applicability of the analysis and to exercise 
and to validate the resulting compter code and (4) a user’s manual for the 
computer code, termed PEPSIS. 
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LIST OF SYMBOLS 


A 

B 

C 

C H 

C P 

D 

D: D 

H 

J 

L 

M 

P 

Pr 

Re 

S 

T 

U 

V 
W 
X 

-► 

V 

h 


Square Matrix of coefficients 

Constant of proportionality 

Constant of proportionality 

Nondimens ional heat transfer coefficient 

Specific heat 

Column vector whose elements are spatial differential 
operators 

Second invariant of the mean flow rate of deformation tensor 
Column vector associated with marching direction terms 
Jacobian 

Linear differential operator 
Mach number 
Static pressure 
Prandtl number 
Reynolds number 

Source term, distance from the leading edge 

Stat ic temperature 

Streamwise velocity component 

Transverse velocity component 

Spanwise velocity component 

Distance from leading edge 

Velocity 

Damping coefficient 

Metric coefficient, enthalpy 

Mixing length 
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q 

w 

X 


y 


a 

Y 

«b 

n 

K 

u 

p 

T 

V 

AX 


i 

J 

Z 

n 

w 

R 

S 


Unit vector normal 

Heat transfer 

Velocity component 

Coordinate direction 

Distance from a surface 

Nondimensional distance from a surface 

GREEK SYMBOLS 


Yaw angle 

Ratio of specific heats 
Boundary layer thickness 
Blasius similarity parameter 
von Karman constant 
Viscosity 
Density 

Viscous stress tensor 
Nabla operator 

Marching direction step size 

Subscripts 

£thd direction 

direction 

Laminar 

Normal direction 

Wall 

Reduced 

Sonic line 
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Subscripts (Continued) 


Turbuent, tangential direction 
Stagnat ion cond it ion 
Streamwise direction 
Cross plane direction 
Cross plane direction 
Free stream condition 

Superscripts 

i ^ streamwise station 


Transpose 



II. ANALYSIS 


Governing Equations : 

The fluid dynamic conservation laws of mass, momentum and energy 
respectively can be written in nondimensional operator form as 


V'/jV = 0 


( 1 ) 


V (/>vv) + Vp - = o 


( 2 ) 


and 




2S . o 

Re 


(3) 


This form of the governing equations, often referred to as the full Navier- 
Stokes equations, requires several auxilliary relationships and models before 
these equations can be solved. In this study, the stagnation enthalpy, h Q , 

is related to the static temperature, T, and the velocity, V, through the 
relationship (assuming constant specific heat) 


h o s C f> T + 


V-V; 


(4) 


while the temperature, T, pressure, P, and density, p, are related by means 
of the calorically perfect gas equation of stae 


y-l 

P = Vc p/ >T 


(5) 


The stress tensor, T, is modelled by the relationship 


r = /i (Vv + VV T ) - -y /xV-V 


( 6 ) 


where the superscript T refers to the transpose of the tensor. The 
components of the velocity vector, V, are interpreted as the mass weighted 
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mean velocity components and p, P and T are the ensemble-averaged density 
pressure and temperature (Ref, 13). Hence, these equations can be applied to 
both laminar and tubulent flows if the effective viscosity, U, is interpreted 
as the sum of the laminar and turbulent, Uj- , viscosities, i.e., 

(7) 

It is assumed that the laminar viscosity can be computed from Sutherland's 
law, and that the laminar and turbulent Prandtl numbers, Pr £ and Pr t are 
constant. For this study, an algebraic mixing length turbulence model of the 
form 


v- T 


Re/oiJy^D 


( 8 ) 


was used where £ m is the algebraic mixing length and D:D is the second 
invariant of the mean flow rate of deformation tensor (Ref. 14). In this 
study, the mixing length of McDonald and Camarata (Ref. 15) was used. 


i m = 0.09 S b tanh[/c7/(0.09S b )]2> 


(9) 


Where 6^ is the local boundary layer thickness, < is the von Karman 

constant, y is the distance to the nearest wall, and 2 ) is the sublayer 
damping term of van Driest (Ref. 16). 

To obtain what is often referred to as the reduced form or the 
'parabolized' form of the Navier— Stokes equations involves approximation of 
the diffusion terms (both stress and Fourier heat conduction) of Eqs . (2), 
(3) and (6). This approximation neglects all derivatives of the stress 
tensor and the Fourier heat conduction terms in a selected 'marching* or 
'streamwise' direction. In addition, all streamwise derivatives of the 
velocity components of the stress tensor are neglected. For example, in a 
general orthogonal coordinate system the principle and shear stress compo- 
nents can be expressed respectively as 



j=. h i h J 


_dhi 

dxi 


_ _L JL y .A. ( -A. w . \ 


j* 


(10) 
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and 


h. A. .Hi _A_ (Ai\ 

T ij P hj dX|\hj/ hj 8xj\h./ 


( 11 ) 


where 


J = (12) 

1 refers to the streamwise direction and 2 and 3 refer to the cross plane 
directions, refers to the velocity component in the i th direction and 
h£ refers to the metric in the i th direction. The approximation neglects 
all direction 1 derivatives of velocity components in Eqs (10) and (11). 
Thus, for example, and are approximated by 




a Ah 4. A± 2. f± \.J_/ A Vl, t , % 

\ h,h 2 dx z h 3 h, dxj 3 j[^x 2 ( h 2 2 ) a *3 ( h 3 W3 )J (1 } 


and 


r „ A a /*i\ 
r -3 = ^ ^r 3 W 


(14) 


Application of the approximations needed to obtain the reduced form of the 
Navier-Stokes equations in other coordinate systems is straightforward. 
Hence, in general, the reduced forms of Eqs. (2) and (3) can be recast as 


V • (/ovv) + VP - 



0 


(15) 


and 




(16) 


where the subscript R refers to the approximated or reduced form of the noted 
term . 
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The reduced form of the Navier-Stokes equations, Eqs . (1), (15) and (16) 
is the starting point for the. discuss ion of the governing equations to be 
used for this study. The intent is to demonstrate that this set of equations 
is not well posed for solution by spatial forward marching when applied to 
the class of problems considered in this study, i.e., supersonic flow with 
embedded subsonic boundary layer regions. Although it does not appear that a 
rigorous analysis has been obtained for the compressible reduced form of the 
Navier-Stokes equations, model sets of equations have been investigated and 
the resuts can be used to give indications of the nature of the compressible 
reduced form of the Navier-Stokes equations. The Euler equations for 
compressible flow are one such relevant model system and are of interest here 
since it is desired to have a stable integration scheme for this system which 
we can reasonably expect to encounter in those high Reynolds number essen- 
tially inviscid regions of the flows considered. It is well known that all 
characteristics of the Euler equations are real for supersonic flow, and thus 
it is inferred that these equations are well posed for solution by spatial 
forward marching (e.g., Ref. 2). Apparently, the supersonic reduced form of 
the Navier-Stokes equations are also well posed for solution by spatial 
forward marching because entirely supersonic flow solutions have been 
obtained using marching techniques (e.g. Ref. 2). In view of the presence of 
imaginary characteristics associated with subsonic flows, it is inferred that 
both the compressible Euler equations (Ref. -17) and the incompressible 
reduced Navier-Stokes equations are unsuitable for solution by spatial 
forward marching (Ref. 18). The fact that these two sets of equations are 
ill posed as initial value problems leads one to suspect that the reduced 
Navier-Stokes equations for mixed supersonic-subsonic flows are also not well 
posed for solution by spatial forward marching. Examination of the 
characteristics analysis for the incompressible reduced Navier-Stokes 
equations shows that the imaginary roots can be affected by the streamwise 
pressure gradient term (Ref. 6). If this streamwise pressure gradient term 
is either a priori specified or neglected, the characteristic equations yield 
only real roots for subsonic flows. Thus, to create a well posed set of 
governing equations suitable for solution by spatial forward marching, much 
effort has concentrated on approximation or modification of this term in the 
reduced form of the Navier-Stokes equations. Efforts to create new sets of 
equations which may be solved by spatial marching which either approximate 
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the reduced Nav ier-Stokes equations or permit stable iterations which, upon 
convergence, represent numerical solutions of the reduced Navier-Stokes 
equations are reported in Refs. 4-8. In view of the present interest in 
flows with strong pressure gradients, equation systems which contain 
approximations to the pressure gradient, in particular the streamwise 
pressure gradient, are viewed with concern. Equation systems which require 
global iteration, i.e. repeated streamwise sweeps through the entire flow, 
yet treat the pressure gradient terms without approximation upon convergence, 
are of course much more preferrable for this class of problems. The present 
desire is to make few approximations to the pressure gradient terms yet 
achieve the computational efficiency of a noniterative forward marching 
algorithm. 

A second set of model equations was investigated by Garvine (Ref. 3). 

In this case, the reduced form of the Navier-Stokes equations is approximated 
by the spatially hyperbolic Euler equations in the inviscid region, and the 
spatially parabolic boundary layer equations in the viscous region. The 
behavior of the solution of this combined set of equations, both of which are 
separately well posed for solution by spatial forward marching in their given 
flow regimes, was shown to be unstable because of the existence of an 
exponentially diverging term in the solution. The unstable behavior was 
caused by the interaction process at the boundary between the two sets of 
equations. A physical interpretation of the unstable interaction process can 
be given for the case of an adverse pressure gradient from a supersonic 
region being impressed on a subsonic region (e.g. a shock wave-boundary layer 
interaction). The adverse pressure gradient causes the subsonic layer to 
increase in thickness. The growing subsonic layer in turn causes the 
supersonic flow to be displaced causing a further increase in the magnitude 
of the adverse pressure gradient. This process is obviously unstable as 
there is no restraining mechanism present at the boundary between the 
subsonic and supersonic regions. The cause of this phenomenon, which is 
sometimes referred to as a departure or branching behavior , has been 
encountered in numerous studies of supersonic interacting boundary layers. 
Note that branching is an entirely different phenomenon than the growing 
modes which cause the subsonic reduced Navier-Stokes equations to be ill 
posed for solution by spatial forward marching. The subsonic reduced 
Navier-Stokes equations possess imaginary characteristics and consequently 
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are ill posed for forward marching, while the two sets of equations analyzed 
by Garvine are individually well posed for forward marching. Rather, it is 
the interaction at the boundary between two sets of equations that causes 
their growing modes. 

The approach taken in this investigation is to find further 
approximations which when utilized within the reduced form of the 
Nav ier-S tokes equations will produce a set of governing equations which are 
well posed for solution by spatial forward marching. It is realized that 
such further approximations to the Navier-Stokes equations will almost 
certainly introduce further limitations on the domain of accurate physical 
represenation of the flow. However, it is believed that approximations can 
be made that will leave the essential physical process of interest intact for 
a wide range of practical cases. It is further believed, that the increase 
in computational efficiency which will result from using a noniterative 
spatial forward marching technique, when compared to current (iterative) 
techniques for solving the full Navier-Stokes equations, justifies the use of 
the approximations . 

The previous discussion of the character of the reduced form of the 
Navier-Stokes equations in mixed supersonic-subsonic flow gives little 
guidance for the choice of further modifications or approximations that will 
yield a set of well posed governing equations. However, three important 
points were made in that discussion: (1) numerical experience indicates that 

the compressible reduced equations are well posed for solution by spatial 
forward marching if the flow is entirely supersonic, (2) the reduced 
equations are known to be ill posed for solution by spatial forward marching 
if the flow is subsonic and (3) the use of different sets of governing 
equations in the supersonic and subsonic flow regions may still be ill posed 
for solution by spatial forward marching even though each set of equations is 
by itself well posed in the region in which it is applied. The problem here 
arises because of the unstable interaction occurring at the boundary between 
the regions. 

In this investigation, the strategy taken is to divide the flow into 
supersonic and subsonic flow regions and to utilize different approximations, 
resulting in different sets of governing equations, in each region. This 
aspect of the approach is not unlike that utilized by other investigators 
(Refs. 5-12), however, the approximations used in obtaining the governing 
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equations in the subsonic region are different than those previously used, 
and hence the interaction of the two sets of equations at the 

subsonic-supersonic boundary is also different. In this study, the technique 
is to utilize the reduced form of the Navier-S tokes equations in the 
supersonic region(s) of the flow, Eqs . (1), (15) and (16), and what can be 
considered to be a model set of equations in the subsonic region(s) of the 
flow. The model set of equations used in the subsonic region(s) is obtained 
by starting with the reduced form of the Navier-Stokes equations and making 
appropriate physical approximations in this region to obtain a new set of 
governing equations such that the coupled system of the inner subsonic flow 
and the outer supersonic flow are stable when solved as an initial value 
problem in space. 

For the problems of interest in this study, high Reynolds number 

supersonic flow with embedded subsonic regions in internal flow devices, the 

boundary layer thickness, 6, will in many cases be small with respect to a 

characteristic vertical dimension of the device. In this investigation, the 

less restrictive assumption is made that the thickness of the subsonic 

portion of the boundary layer is small with respect to the vertical 

dimension. In this subsonic portion of the boundary layer, the usual 

boundary layer approximations for high Reynolds number flow are certainly 

valid. (Note, for instance, that at M * 2 a turbulent flat plate boundary 

layer is supersonic within the viscous sublayer which typically has a 

nondiraensional y + value on the order of 10. Thus, except at very low 

Reynolds numbers, the sonic point is at least one order of magnitude or more 

smaller than the boundary layer thickness). As a result, an order of 

magnitude analysis of the terras in the reduced form of the Navier-Stokes 

equations allows the convection and the diffusion terras to be neglected in 

the subsonic normal (to the wall) momentum equation. This equation can then 

be expressed as a balance between the normal pressure gradient and the 

2 

centrifugal (curvature) forces, 3p/ 3n = U /R, in 2-D streamline coordinates, 
where R the streamline radius of curvature, U the streamwise velocity and n 
the normal to the streamline. In general orthogonal coordinates, Xi, X 2 and 
X 3 with corresponding metric coefficients h\, h 2 and [13 and velocity 
components wi, W 2 and W 3 this equation is expressed as 
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( 17 ) 





Hi'Af.A-, ini'! 

1 1 h , \ 1 ax, * ax , ) 


where X n and X t respectively refer to the appropriate cross-sectional 
direction normal to and tangential to wall (n and t have values of 2 or 3; 
direction 1 is the nominally streamwise direction). It is further possible 
to integrate the continuity equation from the wall to an arbitrary point in 
the subsonic portion of the boundary. This yields in general othogonal 
coord inates . 


h,h T />w n 


= -/ [-^ ( h 2 t V >W l) + 



where again the subscripts n and T refers to the cross flow direction normal 
and tangential to the wall, s refers to the evaluation at the arbitrary point 
in the subsonic region and the subscript w referes to the evaluation at the 
wall. For the class of high Reynolds number flows considered, the boundary 
layer thickness is assumed to be small and as noted earlier the subsonic 
portion of the supersonic turbulent boudnary layer, X s , is usually at least 
an order of magnitude smaller than the boundary layer thickness. Restricting 
our attention to flows where the subsonic region is sufficiently small allows 
the integral in Eq . (18) to be neglected and hence this equation can be 
approximated by 



(19) 


For the case of an impermeable wall Eq . (19) further reduces to 


w n J = o (20) 

s 

In summary then, Eqs . (17) and (18), the streamwise and tangential components 
of the vector Eq . (13), and Eq . (16), the energy equation, constitute the 
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model set of governing equations utilized in the embedded subsonic regions. 

In two-dimensional and axisyrametr ic flows the terms normal and tangential 
directions to the wall are unambiguous and are defined as the X£~ and X 3 - 
directons respectively. For the case where the cross plane is nonaxisym- 
metric (e.g. a rectangular cross section) ambiguity is avoided by referring 
to the normal and tangential direction relative to the nearest wall in the 
subsonic region of the flow. Far from the corners little concern arises from 
use of the nearest wall approach. The corner region is treated by defining a 
corner bisector in the subsonic regions and thus allowing the definitions of 
normal and tangential to change when crossing this line. 

There are several important features of the subsonic model set of 
governing equations. First, no approximation was made to the strearawise 
pressure gradient term (or any other term in the streamwise momentum 
equation). Hence, the full effect of this term will be felt in the subsonic 
portions of the flow. In addition, , the reduced form of the tangential momen- 
tum equation (i.e. the tangential component of Eq . (15)) is unmodified in the 
subsonic region. This allows the effect of the tangential pressure gradient 
to be felt in the subsonic regions as, for instance, would physically occur 
in the case of a glancing shock wave-boundary layer interaction. The assump- 
tion needed to modify the normal momentum and continuity equations in the 
subsonic regions is the relatively unrestrictive condition that the subsonic 
layer is thin relative to the characteristic' transverse dimension of the flow 
device. For the case of an impermeable wall, this leads to the condition 
that within the viscous subsonic layer the normal velocity component is 
negligible, Eq . (20). Since the boundary layer approximation already assumes 
that the normal velocity is small, this condition can be considered to be a 
further approximation, to be applied only in the thin subsonic portion of the 
boundary layer. The importance of the specification of the normal velocity 
is that a mechanism has now been established to prevent the growing mode 
caused by the interaction between the subsonic and supersonic layers, i.e., 
the branching phenomenon. In summary, the new set of governing equations 
consisting of the reduced form of the Navier-Stokes equations in the 
supersonic portion of the flow and the model set of equations in the subsonic 
regions of the flow has, on the basis of numerical experimentation (to be 
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presented herein) been found to be well posed for solution by spatial forward 
marching for a wide range of practical problems. 

Initial and Boundary Conditions : 

To uniquely define the problem of interest, it is necessary to specify 
both initial and boundary conditions. For a spatial forward marching pro- 
cedure, the initial conditions refer to the set of conditiions that must be 
specified at the initial marching station Boundary conditions must be set on 
the boundaries of the cross-sectional marching plane. For the calcualtion of 
internal flows, two types of initial conditions were utilized in this study. 
The first, which is primarily used for flows into devices which have sharp 
leading edges, sets the initial conditions as the free stream conditions. 
Analysis of the characteristics of the supersonic Euler equations shows that 
there are five characteristics entering the upstream boundary of the computa- 
tional domain. Hence, five conditions must be set on this boundary. In this 
study those conditions are chosen as the three velocity components, the pres- 
sure and the temperature. Usually, but not necessarily, the conditions are 
chosen to be uniform everywhere in the initial plane. The second type of 
initial condition is primarily used for cases where information exists at an 
initial plane such that a reasonable approximation to a complete set of ini- 
tial data can be constructed. In its most pure form, this would be an ini- 
tial plane where experimental data were available such that all the initial 
conditions were known. Usually a limited amount of information is available 
where, for instance, free stream conditions, a boundary layer thickness, and 
a skin friction coefficient might be known. In this case a theoretical 
boundary layer profile of the pertinent variables (velocity components, tem- 
perature, pressure, etc.) can often be derived and matched with the free 
stream portion of the flow. It is to be emphasized that the initial 
conditions must in some sense be consistent with the governing equations. In 
supersonic flow computations inconsistencies, perturbations, etc. can persist 
far downstream. 

As used in this investigation, the boundary conditions utilized on the 
bondaries of the cross-sectional plane can be divided into three categories: 
(1) wall conditions, (2) symmetry conditions and (3) external flow condi- 
tions. Analysis of the charcter ist ics of the boundary layer equations shows 
that four conditions must be specified on walls. For this study, the no-slip 
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conditions are used for the strearawise and tangential cross plane velocity 
components, i.e., 

w, = 0 (21) 


and 



( 22 ) 


where again 
subscript T 
cross plane 
normal mass 


the subscript 1 refers to the streamwise direction and the 
refers to the cross plane tangential velocity direction. For the 
normal velocity component either the normal velocity or the 
flux are specified, i.e., 


w 


n 


(23) 


p w 


n 


P n w n 


(24) 


where the subscript w referes to the specified wall value. The fourth 
condition used, the thermal condition, is either to specify an adiabatic wall 
or to specify the wall temperature (a cold or hot wall). The conditions can 
be specified respectively as 


or 


VVT =0 


(25) 


T 


(26) 


where in this case n w represents the unit vector normal to the wall. In 
addition, a fifth condition, not required by the characteristic analysis, is 
used for convenience to close the set of equations. The need for this fifth 
condition could be removed by the use of one-sided differencing or by 
applying one of the governing equations at the wall. In this study, the 
second method was used and the boundary layer approximation to the normal 
momentum equation was applied at the wall. This can be expressed as 


fi w 'VP = 0 


(27) 


Studies have indicated that there is little difference between using this 
equation and the full normal momentum equation (e.g. Ref. 19). 
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The symmetry conditions are meant to be applied on a plane or axis of 
symmetry. The velocity conditions require that the cross plane velocity 
velocity conditions require that the cross plane velocity component normal to 
the axis or plane of symmetry equals zero, i.e., 

n 6 • V = 0 (28) 

where n 8 is the unit vector normal to the axis or plane of symmetry and 
that the first derivatives of the remaining two velocity components equal 
zero. Two other conditions must be set on the axis or plane of symmetry. 
Usually the symmetry conditions on pressure and temperature are used, viz. 

TT $ * VP = 0 (29) 

and 

lf s ‘ VT = 0 (30) 

The final category of boundary conditions used in this investigation 
are those on external surfaces, specifically on the boundary upstream of the 
cowl surface of a supersonic inlet (Fig. 2). In this case a shock wave is 
generated by the ramp and passes out of the computational domain upstream of 
the cowl lip. Upstream of the point where the shock wave passes out of the 
computational domain, the free stream conditions are appropriate as boundary 
conditions and downstream of this point the post shock (Rankine-Hugoniot) 
relationships are valid. Two techniques are commonly used to define shock 
waves, shock fitting and shock capturing. The shock fitting technique 
recognizes the failure of Taylor series expansion through the discontinuity 
so first locates the position of the shock wave, and then enforces the 
Rankine-Hugoniot conditions across the wave. This can occur either on the 
boundary or an interior portion of the computational domain. This technique 
has been used in Ref. 5 to locate the bow shock for external flow cases. 

For internal flow cases, the shock structure can in many cases become very 
complex, and attenuated following the 'boundary layer' interactions. 

Thus, the logic needed to locate the shock waves (especially in three 
dimensions) can become very complex. In addition, in the case of shock 


22 


wave-boundary layer interactions, the shock wave at some point in the 
interaction process ceases being a shock wave and in this region the shock 
fitting procedure becomes very unclear. The shock capturing technique used 
in this study allows the shock waves to be formed as a consequence of the 
solution of the governing equations. Although the presence of the shock wave 
violates the Taylor series representation of the solution used to construct 
the numerical derivatives, it is very convenient. Its use in the present 
study is only justified a posteriori by virtue of the adequacy of the results 
for the problems considered. Returning to the boundary condition upstream of 
the cowl, the approach taken here is to find a set of boundary conditions 
that can be applied on this interior flow 'boundary' that will permit the 
exterior region to be neglected yet that will allow all disturbances which 
originate from within the computational domain to pass through this boundary 
without spurious reflection. The technique used in this investigation is 
predicted on the concept that in a simple wave region the flow properties 
remain constant along Mach lines (the presumption here is that the regions 
fore and aft of the shock wave are simple wave regions). Thus, the first 
derivatives of the flow variables in the direction of the Mach angle should 
be small and are here set equal to zero. The technique is termed Mach wave 
extrapolation and yields the boundary conditions 

r7 • VV =0 (31) 

m 

n_ * VP =0 (32) 

m 


and 


■n m ’ VT =0 (33) 

where is the unit vector in the direction of the local Mach angle. This 
technique requires computation of the Mach angle, and has been successfully 
applied to a number of test cases both by the present authors and the authors 
of Refs. 20 and 21. The boundary conditions allow the flow upstream of the 
shock wave to remain undisturbed, and permit the shock wave to pass out of 
the computatonal domain without reflection. The resulting 'free stream' flow 
behind the shock very closely approximates the appropriate theoretical post 
shock values . 
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III. SOLUTION OF THE GOVERNING EQUATIONS 


The govening equations in both the supersonic and the embedded subsonic 
portions of the flow are simultaneously solved by the consistently split 
linearized block implicit (LBI) technique described in detail in Refs. 2 and 
20. This technique can be logically divided into three parts: 

( 1 ) linarization of the governing equations, ( 2 ) discretization of the 
resulting set of linearized equations by finite difference approximation of 
derivative terras and ( 3 ) simultaneous solution of the resultant set of linear 
coupled algebraic equations. Application of the LBI technique to a set of 
govering equations (and boundary conditions) that is well posed for forward 
marching is straightforward. It is presumed that a solution is known at some 
arbitrary i*^ streamwise station and it is desired to march that solution 
to the i + 1 st station, at some distance AX apart (See Fig. 1). Using 
notation similar to that of Ref. 20 at a single grid point, the system of 
governing equations can be written in the following form: 

= D(<ft) + S(<f>) (34) 

dx 

where <j> is the column vector of dependent variables (w^ , W 2 , W 3 , p, 

M, H and S are column vector algebraic functions of <|>, and D is a column 
vector whose elements are the spatial differential operators which generate 

all spatial derivatives appearing in the governing equation associated with 
that element . 

The solution procedure is based on the following implicit marching 
direction difference approximations of Eq. (34) 

L,i+I ui . , 

H " — = p ** 1 + S 1 * (35) 

Ax 

where, for example, denotes H (<j>^ + *). A local spatial linearization 

(Taylor series expansion about <\> x ) of requisite formal accuracy is 
introduced, and this serves to define a linear differential operator L such 
that 


D , + l = D' 


+ lV +i -V) + 


0(A 


x 2 ) 


(36) 
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Similarly, 


H ' +l ‘HM-f£)V +l - .#>'>+ <W> 

Eqs . (36) through (38) are inserted into Eq. (35) to obtain the following 
system which is linear in <|) 1+ 1 

(A-AxL l )(<£ i + , -<£ i ) * Ax(D' +S') (39) 

and which is termed the linearized block implicit (LBI) scheme. Here A 
denotes a square matrix defined by 

A ■($'-*(§)' 

Eq. (40) is 0 (AX) accuracy. 

It is well known that the finite difference analogue of the governing 
equation system may have an associated stability restriction (Ref. 23). For 
simple equations, the stability criterion can often be analytically obtained, 
for instance using the Fourier technique of von Neumann (Ref. 23). For 
complex systems of equations (including boundary conditions) it is often 
impossible to derive a closed form criterion which can be easily 
interpreted. In this case, the stability bounds (if they exist) may be 
determined by numerical experimentation. For this investigation, second 
order central differences have been used throughout except for the streamwise 
derivatives (which are represented by first order backward differences 
although other choices are clearly permissible). These differences are used 
with the previously described linearized block implicit scheme and a 
numerical analogue of the governing equation system constructed. As far as 
can be determined by experimental investigation, there is no stability 
restriction associated with the resulting scheme. The tests used to 
substantiate these remarks will be discussed later. 

To obtain an efficient algorithm, the linearized system, Eq. (39) is 
split using ADI techniques. To obtain the split scheme, the multidimensional 
oeprator, L, is rewritten as the sum of two 'one-dimensional' suboperators 


(37) 

(38) 
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L£ (i ■ 2,3) each of which contains all terms having derivatives with 
respect to the i^-cross plane coordinate. The split form of Eq. (39) can 
be derived either as in Ref. 20 by following the procedure described by 
Douglas and Gunn (Ref. 24) in their generalization and unification of scalar 
ADI schemes, or using the approximate factorization as in Ref. 26. For the 
present system of equations, the split algorithm is given by 


(A-AxlJh^*-^)* AxtoUs 1 ) 

(41) 

(A -AxLj,)(<£ i + '-<£') =A 

(42) 


where <j>* is the consistent intermediate solution (Ref. 22). If spatial 
derivatives appearing in and D are replaced by the difference formulae, 
as indicated previously, then each step in Eqs. (41) and (42) can be solved 
by a block tridiagonal elimination. 

Combining Eqs. (41) and (42) gives 

(A-AxL|)A" l (A-AxL[)(^> l+l -<#> , ) = Ax(D i +S l ) (43) 

which approximates the unsplit scheme, Eq. (39) to O(AX^). Since the 
intermediate step is also a consistent approxiation to Eq. (39), physical 
boundary conditions can be used for <{>* (Refs. 22, 26). 


26 


IV. TEST CASES 


The preceding analysis was incorporated into a very general computer 
program with the acronym PEPSIS. To validate the ability of this computer 
program to accurately predict flows which are suited for solution by spatial 
forward marching , a series of test cases were run for which there existed 
either experimental data, an analytical solution or a numerical solution 
obtained by another computer program. Calculations were made for both two 
and three dimensional cases in cartesian, general orthogonal, axisymmetr ic 
and nonorthogonal coordinate systems. Both laminar and turbulent test cases 
were considered. In addition, options within the code were constructed to 
solve the conventional two dimensional boundary layer equations, and to solve 
the system proposed by Rudraan and Rubin (Ref. 4) and Lubard and Helliwell 
(Ref. 5) so that results obtained from the present analysis (the reduced form 
of the Navier-Stokes equations) could be compared with other proposals and 
conventional boundary layer calculations. The boundary layer option solves 
the streamwise momentum equation (with specified streamwise pressure 
gradient) and the continuity and energy equations. The numerical solution 
procedure for all these options is the same, appropriately reduced to reflect 
the different sets of governing equations. The only significant numerical 
difference is that in the boundary layer option the continuity equation is 
solved by a trapezoidal integration technique to avoid the need to specify a 
vertical velocity condition at the outer boundary. The boundary conditions 
used in the boundary layer option were the no slip conditions and appropriate 
thermal condition at a wall and specified conditions at the outer boundary. 

CASE I - Incompressible Laminar Flat Plate Boundary Layer 

Initially a low Mach number virtually incompressible zero streamwise 

pressure gradient laminar flat plate case was run with the boundary layer 

analysis to determine how that version of the computer code would reproduce 

the Blasius solution. The case calculated was for a free stream Mach number 

5 

of 0.1 and a Reynolds number per unit length of 10 per meter. The wall 
temperature was chosen as the adiabatic wall temperature. The computational 
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domain is as shown in Fig, 3. 100 grid points were nonuniformily distributed 

in the transverse direction with grid points concentrated in the region close 
to the wall. The initial boundary layer profile was generated from a Blasius 
solution at an axial location of X/L = 2.0 from the leading edge of the flat 
plate. (L was chosen to equal 1 meter). The initial boundary layer (point 
where u/ue = 0.999) was contained within the first 65 grid points 
corresponding to a thickness of 6/L = 0.0267. The initial displacement and 
momentum thickness Reynolds numbers were 297.0 and 771.7 respectively. The 
initial profile was marched downstream 500 steps to a streamwise location of 
X/L = 6.99 at a constant marching step size of AX/L = 0.01 (this corresponds 
to a streamwise marching step size of Ax/ 6 = 0.375 of the initial boundary 
layer thickness). The calculated streamwise velocity profile at X/L * 6.99 
(plotted in terms of the Blasius similarity variable n = Y / pu«/ M«X) 
iscompared with the theoretical Blasius profile in Fig. 4. Agreement with 
the Blasius profile is excellent. In Fig. 5 the calculate^ streamwise 
distribution of skin friction coefficient is compared with the Blasius 
result. Agreement between the two results is again excellent. The ability 
of the boundary layer version of the code to accurately predict the Blasius 
solution is viewed as a prerequisite before more complex cases can be 
attempted. In addition, since the subsonic layer approximations used in the 
more general analysis are similiar to the boundary layer equations, the 
numerical scheme must be able to solve this related set of equations. 

CASE II - Supersonic Laminar Flat Plate Boundary Layer 

Next both the boundary layer option and the more general analysis were 
used to predict the laminar supersonic flow over a flat plate. This case was 
confined to a model fluid having a laminar viscosity proportional to the 
temperature and a unity Prandtl number. For this case, the 

Dorodnitsyn-Howar th similarity solution of Ref. 27 can be used as an initial 
condition and as a means of generating a theoretical downstream solution. 

The case run was for a free stream Mach number of 5.0 the Reynolds number per 
unit length was 10 5 per meter, the wall temperature was specified at 25% of 
the free stream temperature and the reference length, L, was 1.0 meter. 
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One hundred (100) grid points were nonuniformly distributed in the transverse 
direction with grid points concentrated in the region close to the wall. The 
initial profile was specified at a streamwise location of X/L = 2.0 from the 
leading edge. The initial displacement and momentum thickness Reynolds 
numbers were 297.0 and 771.7 respectively. The initial boundary layer 
profile was contained within the first 68 grid points corresponding to a 
boundary layer thickness of S/L = 0.0359. The subsonic portion of the 
initial boundary layer was contained within the first 25 grid points 
corresponding to a subsonic thickness of 6 . 55 % of the boundary layer 
thickness. Both the boundary layer and the more general analysis versions of 
the computer code ^ere utilized to march the initial solution downstream in 
500 equal steps of AX/L = 0.01 (corresponding to a Courant number of 0.81) to 
a streamwise location of X/L = 6.99. The boundary conditions for the 
boundary layer version were the same as for the previous case except that the 
wall temperature was now specified. The more general analysis utilized the 
no slip conditions, specified zero normal pressure gradient and specified 
temperature at the wall and Mach line extrapolation at the outer surface as 
boundary conditions. 

The streamwise velocity profiles (plotted at X/L = 6.99 in terms of the 
Blasius similarity variable h = y /pucV UooX) predicted by both the bound- 
ary layer and the more general analysis are compared with the theoretical 
profile in Fig. 6 . As can be seen, the agreement is excellent as the 
calculations and the theoretical solutions are indistinguishable from one 
another. In Fig. 7 the calculated streamwise distribution of surface skin 
friction coefficient is compared with the theoretical values. As can be seen 
the agreement between the boundary layer and more general analysis and the 
theoretical distribution is excellent. Similar excellent agreement can also 
be observed in the plot of the streamwise distribution of momentum thickness 
Reynolds number, Re 0 (see Fig. 8 ). 

The purpose of the above case was to demonstrate that when the interac- 
tion effects of displacement are negligible the more general analysis yields 
approximately the same results as a boundary layer analysis (in this case the 
boundary layer profile, a local property - the skin friction coefficient and 
an integrated property - the momentum thickness Reynolds number were com- 
pared) . The inference is, of course, that the assumptions utilized to make 
the governing subsonic equations of the more general analysis well posed for 
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solution by spatial forward inarching did not compromise the physics of the 
flow. The case chosen was a high Mach number flow with a highly cooled wall 
where the displacement effects were expected to and evidently did have only a 
small influence on the boundary layer flow. 

CASE III - Supersonic Turbulent Flat Plate Boundary Layer 

The next case considered was the supersonic turbulent flow over a flat 

plate. As with the previous case solutions were obtained with both the 

boundary layer and the more general analysis. The energy equation was 

approximated by assuming constant stagnation enthalpy. Boundary conditions 

were the same as used for the previous case. The method of Maise and 

McDonald (Ref. 28) was applied to the incompressible Musker profile (Ref. 29) 

to obtain both an initial compressible flow condition and as a basis for 

comparison with the calculated downstream results. For this test case, the 

free stream Mach number was chosen as 3.0, the Reynolds number per unit 

length was 10 5 per meter and the reference length, L, was 1.0 meter. An 

initial boundary layer thickness of 6/L = 0.1365, a momentum thickness 

. 3 

Reynolds number of 934.6, and a skin friction coefficient of 2.5 * 10“ 
were assumed. 50 grid points were nonuniformly distributed in the vertical 
direction with the initial boundary layer encompassing 17 grid points. The 
subsonic portion of the initial boundary layer had a thickness of 0.01152L 
and was contained within the first 4 grid points. The corresponding subsonic 
nondimens ional distance y + = Pwy u T/l*w was 9.5. The initial profile 
was located at a value of X/L =1.0 and was marched downstream in 350 unequal 
steps to a downstream location of X/L = 180. The initial step size was AX/L 
= 0.15 corresponding to a Courant number of 2.88. 

Results from the test case are presented in Figs. 9-12. In Figs. 9 and 
10 the skin friction coefficient vs. the momentum thickness Reynolds number 
results generated by both versions of the computer code are compared with 
results from the transformed profile of Musker. Except for some relatively 
minor deviations near the initial station the agreement is good for both 
cases. It is believed that the minor deviations in this region are due to 
the numerical method adjusting to the given initial profile (which for 
convenience assumed zero initial transverse velocity). The transformed 
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generalized velocity defect predicted by both versions of the code was 
compared with those given by the analytic Musker profile in Figs. 11 and 12 
respectively. As can by seen from these figures the theoretical and 
calculated agreement is good. 

CASE IV - Hypersonic Laminar Corner Flow 

The fourth case considered was that of the hypersonic laminar strong 
interaction flow in a 90° corner formed by two sharp flat plates aligned with 
the free stream. A schematic of the flow system is shown in Fig. 13 along 
with the prescribed coordinate system. A viscous layer starting near the 
leading edge forms the continuum merged layer; the strong interaction regime 
appears downstream with a discrete boundary layer, inviscid region and shock 
wave structures. In the corner region the two layers which form on each of 
the plates merge together and it is this region in particular that is 
examined here. The computational study was conducted at a free stream Mach 
number of 11.2 and a Reynolds number of 5.9 x 10 5 per meter. The reference 
length, L, was chosen as the height and width of the computational domain, 
0.134 meters. This case was experimentally studied by Cresci (Ref. 30). The 
free stream and wall temperature were 361°K and 305.55°K respectively. To 
determine the distribution of cross plane and streamwise grid points 
necessary to adequately resolve the physics pf this case, the two dimensional 
analog of this case was first run, i.e., hypersonic laminar flow over a flat 
plate. Using criteria determined from running the two dimensional case, it 
was decided that the cross plane would require a 50 x 50 mesh of grid 
points. Grid points were packed in the vicinity of the walls and the shock 
region (see Fig. 14). The boundary layer option is inappropriate for this 
and subsequent cases and so was not run. The more general analysis was 
forward marched 120 streamwise steps corresponding to a streamwise location 
of X/L = 1.316. For the first 100 steps (corresponding to a Courant number 
of 0.286) the step size, AX, was kept constant at AX/L = 0.01; thereafter the 
step size was allowed to increase by 5% per step. In the calculation, 
uniform free stream conditions were used as initial conditions. After 
marching the free stream conditions for two streamwise steps, the flow 
encountered the leading edge of the corner and the flow was allowed to 
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naturally develop. In this approach the leading edge singularity is ignored 
and smeared over by the computational scheme. The governing equations used 
for this case were the three momentum equations, the continuity equation and 
the energy equation. The no-slip conditions, the zero pressure gradient 
condition and specified temperature were used as the wall boundary 
conditions. On all other boundary surfaces the symmetry conditions were 
imposed . 

Fig. 15 presents the comparison of the calculated and experimental wall 
pressure distribution normal to either of the two flat plates (i.e. in the Z 
direction) at streamwise location X/L = .990 corresponding to an interaction 
parameter, x = 5.1 (X“ /Roo^^, R«> = P» UooX/ U» where C 

is the constant of proportionality between viscosity and temperature and X is 
the distance from the leading edge). It can be seen that x a 1/X^/^, 

Ref. 31. Considering the uncertainties in the experimental data, the 
calculated wall pressure agrees very well with the measurements both in the 
location of the peak pressure and the general form of the pressure 
distribution. Fig. 16 compares the Stanton number C H where 
C H = q w / [ PooUco( h - h w ] ) at streamwise location X/L = or X = In 

general the agreement is good, with the analysis showing a slightly thinner 
peak heating region. Fig. 17 shows the comparison between computed and 
measured skin friction coefficient at X/L = 0.990 or X ” 5.17. The apparent 
discrepancy between the data and the calculation is large, although the 
general form of the curves are similar. Since it was not clear why such a 
large discrepancy occurred, the original source of the data (Ref. 30 and 32) 
were reviewed carefully to determine possible sources for such disagreement. 
In Ref. 30 it is reported that the skin friction coefficient was estimated by 
using the gradient of axial velocity normal to the wall which was calculated 
from the measurements of the total temperature and Mach number distributions 
in the corner regions. Furthermore, it was found that the nearest measuring 
station was approximately 0.1 cm from the wall. Thus, for the purpose of 
comparison with the data, it was felt to be reasonable to calculate the 
predicted skin friction coefficient in a like manner by using the numerically 
calculated velocity gradient normal to the wall with the first grid point 0.1 
cm off the wall. As is shown in Fig. 17 (dotted line) this calculation 
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produced a skin friction coefficient distribution that is in much better 
agreement with the experimental data. Fig. 18 shows the streamwise 
development of the stagnation pressure isobars. The development of the 
merged region, separation of the shock wave from the viscous region and the 
complex corner flow structure are easily identified as the flow proceeds 
downstream. 

CASE V - Three-Dimensional Glancing Shock Wave - Boundary Layer Interaction 

A typical flow phenomenon occurring in internal supersonic flow is the 
interaction of a sidewall boundary layer with a cowl generated glancing shock 
wave. This flow configuration gives rise to a strong interaction region in 
the corner between the sidewall and the cowl resulting in the formation of a 
corner vortex as the flow proceeds downstream. A well documented extensive 
experimental investigation of the phenomenon has been performed by Oskam, 

Vas and Bogdonoff and is reported in Ref. 33 and Ref. 34. Fig. 19 sche- 
matically depicts the flow. A supersonic turbulent boundary layer is pro- 
duced on the walls of the test section. A shock generator in the form of a 
sharp edged plate is mounted vertically between the tunnel floor and ceiling 
and turned to some desired angle, <f>, to the incoming flow. The glancing 
shock formed by the generator then interacts with the boundary layer formed 
along the floor of the test section. For the case under consideration in 
this study the plate was inclined at <(> = 10° to the free stream which had a 
Mach number of 2.94. The pre-interaction boundary layer thickness was 1.40cm 
and the Reynolds number based on that thickness was 9.68 x 10 5 . For this 
case experimental data were obtained at three stations (see Fig. 20). The 
mean flow data taken were static pressure (cone-cylinder proble) , stagnation 
pressure (cobra probe), total temperature ( thermocouple probe) and yaw angle 
(cobra probe). The yaw angle, a, in the context of the coordinate system of 
Fig. 20 is defined as 


a 


= tan 



(44) 
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where u is the streamwise velocity and v is the velocity component 
perpendicular to the wind tunnel side wall. 

The cross-sectional area of the test section in 8 x 8 inches, and hence 
the reference length, L, was chosen as 8 inches (20.32 cm). The 
computational domain chosen for this calculation consists of the region 
starting at X = 2.98 cm, i.e., X/L = 0.147 upstream of the leading edge of 
the shock generator and proceeding to a distance of X = 25.14 cm, i.e., 

X/L = 1.237 downstream of the leading edge. Because of the vertical 
symmetry, the computation only had to be made in the lower half of the test 
section. The s panwise domain was bounded by the shock generator on one 
surface, and a free stream boundary located far enough away such that the 
shock wave will not exit through this surface for the streamwise extent of 
this computation. The coordinate system used for this calculation was 
generated by using a Schwarz-Christof fel transformation technique of Anderson 
(Ref. 35) to generate a set of conformal coordinates. The initial conditions 
were calculated by assuming a mixture of a boundary layer profile on the 
floor of the test section and free stream conditions elsewhere. The initial 
boundary layer velocity profile was calculated by the method of Maise and 
McDonald (Ref. 28) with a transverse velocity set to zero. The 
pre- interaction boundary layer thickness of 1.40 cm ( 6/L * 0.0689), a skin 
friction coefficient of 1.2 x 10” , a momentum thickness Reynolds number of 
3 x 10”\ and a wall temperature of 297° were used to calculate the 
profile. The initial enthalpy and temperature profiles were calculated by 
use of the modified Crocco profile (Ref. 36). Boundary conditions on the 
shock generator surface, and the floor of the test section are the no-slip 
conditions for the three velocity components, specified temperature for the 
terraal condition (279°K for the test section floor and 236°K for the shock 
generation surface), and imposition of the normal pressure gradient equals 
zero condition. The two above surface temperatures are nominally average 
values of these parameters during a run. Due to the nature of test 
facility, the free stream stagnation temperature decreases on the order of 
55 °K during a run. This causes the surface temperature to also vary; 
however, since the shock generator has a low heat capacity, its temperature 
decreases more rapidly than the test section floor temperature, and thus the 
above temperatures are the mean temperatures during the fun. This 
time-dependent nature of the experiment will undoubtedly lead to some 
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(undetermined) error in assessing the values of the experimental data. The 
boundary conditions on the bounding surface upstream of the leading edge of 
the shock generator as well as the two other boundary surfaces (the plane of 
symmetry and the surface opposite the shock generator) utilize symmetry 
boundary conditions. The mixing length model is based on minimum distance to 
the nearest wall, in this case the distance from a surface is taken to be the 
minimum of the distance to either the shock generator surface or the test 
section floor. The laminar and turbulent Prandtl numbers used for this 
calculation were 0.74 and 0.90, respectively. The initial solution was 
marched downstream 400 stations at a computational step size that varied from 
a minimum of AX/L = 0.002 in the vicinity of the tip of the shock generator 
to a maximum of AX/L = 0.005 further downstream. At the initial plane, this 
minimum step corresponds to a physical distance of approximately AX/L - 0.02 
or a Courant number of 8.70. The cross plane utilized a 40 x 40 grid point 
strucutre with grid point packing about both the shock generator surface and 
the test section floor. As was the case for the previous calculation, the 
shock wave was generated as part of the solution rather than as part of the 
initial profile. The equations solved in this case were the three orthogonal 
momentum equations, the continuity equation and the energy equation. 

Results in the form of calculated and measured pitot pressure, static 
pressure, total temperature and yaw angle are presented in Figs. 21-29. 

Pitot pressure measurements were obtained at the four measuring stations 
shown in Fig. 20, static pressure and total temperature at two stations and 
yaw angle at three stations. Agreement between the calculations and the data 
is excellent both from a qualitative and quantitative perspective. (Data at 
a value of YG = .25 in is within the wind tunnel floor boundary layer. 

All other data is outside the boundary layer.) The agreement between the 
calculated and experimental static pressure (Figs. 21 and 27) indicates 
proper placement of the shock wave. Results at the streamwise station 
X = 7.60 in vertical distances from the shock generator of 2.75 and 
above (i.e. values of YG (Fig. 20)) there is some minor deviation 
between the calculation and experiment. This is probably due to some 
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smearing of the shock wave as it was in this region that the fewest number of 
transverse grid points were used. The pitot pressure profiles (Figs. 23-26) 
again show overall good agreement between the calculation and data. The 
pitot pressure can be viewed as a composite variable which measures the level 
of the static pressure as well as the boundary layer streamwise velocity 
profile. The near wall data values of pitot pressure expecially in the 
corner region show some deviation from the calculated values. Comparisons 
between data and calculated values of the yaw angle are shown in 
Figs. 27-29. The agreement is excellent both qualitatively and 
quantitatively throughout. Some minor deviations can be seen within the 
boundary layers. However, when one considers that it is usually more 
difficult to accurately calculate the cross flow velocity components than the 
streamwise velocity component, the amount of deviation must be considered 
minimal. The yaw angle distribution is one means of determining both the 
position and strength of the corner vortex which is located in the corner 
region and grows as the flow proceeds downstream. The ability to accurately 
calculate the yaw angle, therefore, implies an ability to accurately 
calculate the strength and location of the vortex. Referring back to the 
statement about the local maximum in the total temperature, it can be seen 
from examining the yaw angle distribution that the local maximum in total 
temperature occurs in the vicinity of the edge of the vortex associated with 
the shock wave. 

CASE VI - Supersonic/Turbulent Flow in a Variable Area Ratio Duct 

The next test case considered was the two-dimensional turbulent flow 
through a channel of varying cross-section. The geometry and computational 
mesh used for this calculation is shown in Fig. 30. For this case, a 
nonor thogonal coordinate system was utilized. The streamwise coordinate is 
obtained by using a Y/6(X) transformation, where 6(X) is the equation of the 
height of the top surface taken normal to the lower, flat surface. The 
equations solved were the transformed Cartesian streamwise and transverse 
momentum equation and the transformed continuity equation. The stagnation 
enthalpy was assumed constant. Boundary conditions on both surfaces are the 
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no-slip conditions and zero normal pressure gradient. The flow conditions 
for this case were a free stream Mach number of 1.9 and a Reynolds number per 
unit length of 4.64 x 10 per meter. The reference length, L, is taken as 
the minimum distance between the upper and lower surfaces, i.e., 0.01018 
meters. Ninety-nine (99) grid points were utilized in the transverse 
direction with packing in the vicinity of the two walls. The marching step 
size was taken at a constant value equal to AX/L * 0.01 (corresponding to a 
Courant number of 1.37). The initial profiles on both surfaces were again 
generated by using the method of Maise and McDonald (Ref. 28). The boundary 
layer thickness on the lower and upper surfaces were 6 ]/L = 0.06074 and 
62 /L = 0.12149, respectively. The corresponding skin friction coefficients 

_3 3 

were 3.94 x 10 and 3.32 x 10 . Corresponding momentum thickness 

Reynolds numbers were 400 and 750, respectively. These values were the same 
as were used in the calculation of Ref. 38. In Ref. 38 the same prediction 
was performed by numerically solving the full Navier-Stokes rather than with 
a reduced form as is done here. 

Figure 31 shows the comparison of the upper surface static pressure 
distribution calculated by both the spatial forward marching analysis and the 
Navier-Stokes code. Initially, the agreement between the two predictions to 
good, but the forward marching procedure predicts a higher peak pressure in 
the compression region and a corresponding greater expansion further 
downstream. The differences of the two methods in this region could be due 
to the neglecting of the strearawise diffusion in the forward marching 
analysis (and hence the neglecting of upstream influence). The difference 
could also be due to the higher level of accuracy obtained by the forward 
marching procedure due to the use of a finer streamwise grid structure than 
was used in the Navier-Stokes analysis. The forward marching run was 
terminated as the shock wave approached the lower wall due to the formation 
of a suddenly expanding subsonic region near the lower wall. The streamwise 
location that this phenomenon occurred was approximately the same location 
that the Navier-Stokes analysis predicted the existance of a Mach stem, and 
hence a corresponding suddenly expanding subsonic region. Although the 
forward marching procedure cannot calculate through such a region, it is 
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significant that at least in this case, the analysis predicted the initial 
formation of such a region. 

VII - Supersonic Turbulent Shock Wave-Boundary Layer Interaction 

A well-documented experimental investigation of a shock wave-boundary 
layer interaction flow has been made by Rose (Ref. 36). A schematic of the 
experimental apparatus is shown in Fig. 32. A conical shock wave is 
generated by a 9° half angle cone situated in the center of an axisymmetric 
test section of radius, L - 2.64 cm. The shock wave interacts with a 
turbulent boundary layer on the wall of the test section. The free stream 
Mach number was 3.88 and the Reynolds number based on the pre-interaction 
boundary layer thickness of 0.51 cm was 8.7 x 10 4 . Experimental data were 
obtained in the turbulent boundary layer on the wall of the test section in 
the vicinity of the interaction region. The mean flow data consists of 
measured pitot and total temperature profiles and surface static pressure 
distribution. The cone half angle was chosen to produce a shock strength 
near to that required to produce streamwise separation. 

The computational domain for this calculation consists of the transverse 
region between the cone surface and the test section wall starting at the 
streamwise position X/L = 0.364 upstream of the cone tip and extending 
downstream X/L = 5.053. Because the cone is placed in the center of the 
axisymmetric test section the resulting flow is axisymmetric, and hence it is 
only necessary to solve the usual axially symmetric set of governing 
equations. A conformal coordinate system was generated by means of the 
previously discussed Schwarz-Christof fel transformation technique of Anderson 
(Ref. 35). The initial conditions consist of a turbulent boundary layer on 
the test section wall generated by the method of Maise and McDonald 
(Ref. 28). The boundary layer thickness was 0.51 cm, ( d/L = 0.1923), the 
skin friction coefficient was 1.72 x 10~ , the momentum thickness Reynolds 
was 2000, and the wall temperature was 277. 8°K. The initial enthalpy and 
temperature profile is calculated using the modified Crocco— Busraann profile 
(Ref. 37). Finally, the pressure on the initial plane is assumed constant at 
the test section free stream value. The equations solved were orthogonal 
axisymmetric streamwise and transverse momentum, continuity and energy 
equations. Boundary conditions on the cone surface and the test section wall 
are identical, i.e. no-slip for the velocity components, a specified wall 
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temperature of 277. 8°K and zero normal pressure gradient. On the axis of 
symmetry upstream of the cone tip, symmetry conditions are imposed. The 
laminar and turbulent Prandtl numbers were set to 0.71 and 1.0, 
respectively. The initial solution was marched downstream 800 stations at a 
constant computational step size of AX/L = 0.01 corresponding to a Courant 
number of 0.42. At the initial station this corresponds to a physical step 
size that is 1% of the distance from the tip of the cone to the test section 
wall. For this calculation, 99 transverse grid points were utilized with 
packing about the cone and test section wall. It is to be noted that the 
shock wave generated by the cone is not input as an initial condition. 

Rather, the shock wave is generated as a result of the coordinate system, the 
governing equations and the applied boundary conditions. 

Results in the form of calculated and measured pitot pressure are 
presented in Fig. 33. The experimental data were obtained at equally spaced 
streamwise stations in the vicinity of the interaction. Basically, the first 
three profiles are in the pre-interaction region and the remaining are in the 
interaction and post interaction regions. The incident shock wave can be 
seen in the pitot pressure profile plots (Fig. 33), at the second and third 
data stations. At the sixth and subsequent data stations, the shock wave 
reflects off the boundary layer and proceeds back towards the cone. As can 
be seen from Fig. 33, the calculated and measured values of pitot pressure 
are in substantial agreement. Qualitatively, they agree at all data 
stations. Qualitatively the agreement in the pre-interaction and throughout 
most of the interaction region is excellent. The calculation does predict a 
slightly thicker (on the order of 10%) emerging boundary layer thus resulting 
in some disagreement between the calculated and experimental pitot pressures 
in the outer portions of the boundary layer at the downstream stations. This 
can perhaps be attributed to the use of a constant boundary layer thickness 
in the turbulence model of Eq. (9). To date, no attempt has been made to use 
a varying boundary layer thickness in the turbulence model. Initially, the 
wind tunnel boundary layer had four grid points within the subsonic portion 
of the layer ( correspond ig to 0.004% of the distance between the wall and the 
axis of symmetry). During the interation process, the adverse pressure 
gradient causes the flow to decelerate. At one point the subsonic portion of 
the boundary layer was contained within 12 grid points corresponding to 2.1% 
of the distance between the cone and the wind tunnel wall. 
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Numerical simulation of the Rose experiment were also performed with the 
computer code constructed to perform on option the calculation according to 
the method of Rudman and Rubin (Ref. 4) and Lubard and Helliwell (Ref. 5). 

It will be recalled that Rudman and Rubin neglected the streamwise pressure 
gradient in their procedure while the Lubard and Helliwell technique used 
explicit spatial lagging for the calculation of the streamwise pressure 
gradient term. The Rudman-Rubin approach was not proposed for use in this 
type of flow (it was proposed for hypersonic flow over highly cooled walls 
where the streamwise pressure gradient is negligible with respect to the 
other terras of the momentum and energy equations). Here it is used only to 
give perspective to the role that streamwise pressure gradient term plays in 
the flow of interest. For both techniques the linearization of the governing 
equations, their finite difference approximation and solution of the 
resultant set of linear algebraic equations was achieved by the same method 
as described in the previous section, i.e., the same method as was used to 
solve the governing equations of the present analysis. Although this 
numerical method of solving the governing equations was not used by either 
Rudman and Rubin or Lubard and Helliwell, the governing equation and finite 
difference representation of derivatives, grid point distribution, boundary 
conditions and initial conditions were the same. Hence, the only factor of 
consideration is the governing equations themselves. 

Results for the simulation of the Rudman and Rubin techniques are 
presented in Fig. 34. At this first streamwise data station the calculated 
boundary layer profile agrees well with the experimental data as they should 
since at this point, i.e., in the preinteraction region, the streamwise 
pressure gradient is small. However, at the third streamwise station the 
calculated shock wave location noticeably lags the data, and as the flow 
proceeds downstream the calculated flow bears little resemblance to this 
data. It should be re-emphasized here that the Rudman and Rubin technique 
was not developed for this problem and that the results indicate only that 
this type of approach is not satisfactory for predicting typical flow 
phenomena that occur in internal flow devices. 

Using the method of Lubard and Helliwell no stable calculation could 
be obtained with any marching step size used. The calculations displayed 
the well known symptoms of branching, i.e., a large increase in pressure 
was predicted which in turn causes a strong streamwise recirculation zone 
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to form. Initially this calculation was run with a streamwise marching step 
of 0*01% of the reference length, L (the same as was used for the general 
analysis and the Rudman and Rubin calculation)* Subsequent calculations 
using a marching step as large as 10% of the reference length were attempted, 
but these calculations were also unstable. The step size of 10% of the 
reference length corresponds to approximately twice the stable marching step 
size predicted by Lubard and Helliwell in Ref. 9. This is not surprising as 
the authors infer that their criterion is only approximate. No step sizes 
larger than 10% were used as it was felt that streamwise step sizes of this 
magnitude are surely too large to adequately resolve the interaction 
process* The failure of this method to give a stable solution supports the 
previously stated objection that methods that have a minimum step size 
criterion may not be acceptable for flow situations of interest where large 
streamwise pressure gradients exit. It is to be expected, although it has 
not been demonstrated here, that the other methods which have a similar 
minimum step size criterion; would also not give a stable solution for this 
case . 

VI - Axisymmetric Inlet 

The last two test calculations performed under this effort were for the 
Boeing axisymmetric mixed compression type inlet. Details of this inlet and 
the experimental test conditions are available in Ref. 39. A schematic of 
the inlet is provided in Fig. 35. The test conditions were for a free stream 
Mach number of 3.5 with a Reynolds number (based on the lip diameter) of 
2.8 x 10“ 6 . The reference length L, was chosen as half of the lip diameter 
of 49.723 cm. Calculations were performed at both 0° and 3° angle of 
attack. Surface contours of the centerbody and cowl are listed in Table 1. 

At the design Mach number of 3.5, the centerbody is in the fully retracted 
position providing a capture mass-flow ratio of unity. The axisymmetr ic 
centerbody' s half angle is initially inclined at 10° to the horizontal. The 
Schwarz-Christof fel transformation technique of Anderson (Ref. 35) was used 
to generate a conformal coordinate system for this inlet (see Fig. 36). The 
computational domain consisted of the region between the centerbody and the 
cowl. The upstream limit of the computational domain was chosen to be 
slightly upstream of the leading edge of the centerbody; the downstream 
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extent terminated downstream of the geometric throat. For the 0° angle of 
attack case the flow is axisymmetric and hence only the axisymmetr ic 
equations had -to be solved. In this case, the equations solved were 
streamwise and transverse momentum, continuity and the energy equations. 
Ninety-nine (99) grid points were used in the transverse direction with grid 
packing in the vicinity of the centerbody and the cowl surface. For the 3 
angle of attack case, the flow is three-dimensional. However, a plane of 
symmetry exists and hence the computation only has to be made in the half 
plane. In this case, streamwise, radial and circumferential momentum, 
continuity and energy equations had to be solved. Fifty radial grid points 
were used for the 3° case with grid packing in the vicinity of the centerbody 
and the cowl surfaces. Nineteen equally spaced grid points were used in the 
circumferential direction (corresponding to 10° increments). The boundary 
conditions on both the centerbody and the cowl surface were the no slip 
conditions for the velocity components, zero normal pressure gradient and the 
adiabatic wall condition. The effects of bleed (which was utilized in the 
experiment) were not considered for these two cases. On the free surface 
corresponding to the upstream extension of the cowl surface, Mach line 
extrapolation was utilized. Since the calculation procedure was initiated 
upstream of the leading edge of the centerbody, a boundary condition had to 
be set on this surface (corresponding to the upstream extension of the 
centerbody surface). In this case, symme try .conditions were used for 
streamwise velocity (circumferential velocity for the 3 angle of attack 
case), pressure and temperature. The normal velocity component on the 
upstream extension remained unchanged. Finally, for the 3 angle of attack 
case symmetry conditions were used on the plane of symmetry. The initial 
condition for both cases were the uniform free stream conditions. Thus, the 
shock wave is produced by the calculation procedure and is not input as a 
part of the initial conditions. 

For the 0° angle of attack case the initial conditions were marched 470 
streamwise steps at a constant step size of AX/L — 0.02. Corresponding to a 
Courant number of 3.00. This corresponds to an initial step size equal to 2% 
of the distance from the centerbody to the upper surface (see Fig. 35). The 
computation terminated slightly upstream of the geometric throat (a large 
recirculation zone formed on the cowl surface). The displaced scale Mach 
number profiles are shown in Figs. 37 and 38. Figure 37 shows the plots in 
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physical space while Fig. 38 shows the plots in computational space. 

Although the propagation of the shock waves is somewhat difficult to discern 
in Fig. 37, the shock waves propagation is very distinct in the computational 
space plots of Fig. 38. The centerbody shock passes out in front of the 
cowl. A shock wave forms off the cowl, impinging on the centerbody, 
reflecting and impinging on the cowl surface where a recirculation zone 
forms. Comparison of the centerbody and cowl static pressure distributions 
are presented in Figs. 39 and 40. Considering that the effects of bleed are 
not considered, the agreement with data is good. In Fig. 38, it can be seen 
that the region where the calculation terminated due to a large recirculation 
zone was also (evidently with good reason) a region where the flow was 
subjected to wall bleed. 

The 3 angle of attack case was run mainly as a demonstration case of an 
off-design condition. The centerbody location was specified to be that of 
the 0 angle of attack location. Under this condition, it is to be expected 
that the shock wave will hit inside the cowl on the windward side and fall 
further outside the cowl on the leeward side. However, the strength of the 
shock will be stronger on the windward side. This case was run with a 
variable marching step size. The step size was chosen such that initially a 
step size of Ax/L = 0.02 was used (corresponding to a Courant number of 
1.45). After encountering the centerbody the step size was gradually 
increased to a step size of AX/L = 0.4 and then gradually decreased to a step 
size of AX/L = 0.01 slightly upstream of the cowl. Downstream of the cowl 
entrance region the step size was gradually increased again to a maximum of 
AX/L =0.04. The purpose of varying the step size is to increase resolution 
in regions of large streamwise flow gradients and to decrease resolution of a 
larger step can be taken. This calculation was marched downstream 270 steps 
before the calculation terminated. Termination was due to the generation of 
a recirculation zone on the leeward side of the centerbody. 

Results of the calculating of the off-design 3° angle of attack case are 
presented in Fig. 41-48. No comparison is made with data since data was not 
obtained for the 3° case in the off-design condition. In Figs. 41-43, the 
streamwise pressure distributions on the centerbody and the cowl surfaces are 
shown on the leeward, waterline and windward rays. On the leeward side of 
the cowl, the shock passes in front of the cowl entrance. The turned flow 
then encounters the cowl entrance and forms a shock wave. This can be seen 


43 


by the spike-like rise followed by a decrease in pressure (as the cowl at 
this point becomes concave) followed by an increasing rise in pressure as the 
cowl surface becomes convex again. Similar features can be seen on the 
waterline and windward rays. In this case the shock formed by the centerbody 
impinges on the cowl surface. A similar spike-like behavior followed by 
expansion and compression zones can be observed. The corresponding displaced 
scale Mach number plots are shown in Figs. 44-46. By looking up into the 
inlet, the shock structure can be discerned. Figure 47 shows typical 
secondary flow velocity vectors. The streamwise location is approximately at 
a value of X/L = 2.45. The winward side is on the left; the leeward side is 
on the right. The position of the shock wave is clearly discernable as a 
discontinuity in flow direction. It can be seen that the shockwave is closer 
to the windward surface than to the leeward surface. Mach number contours at 
the same station are shown in Fig. 48. The approximate position of the shock 
wave is shown by the concentration of the Mach number contours. 

All of the above test cases were run on the NASA— LRC IBM 370—3033 

Q 

computer. The CPU run times are 5.43 x 10” sec/grid point for the 
two-dimensional runs 1.43 x 10 -2 sec/grid point for the three-dimensional 
runs. The total run time for a given case scales linearly with respect to 
the number of grid points, thus for instance the Rose case which used 99 
transverse grid points and marched 800 streamwise stations had a run time of 
430 seconds. The difference in run times between two- and three-dimensions 
is due to the three-dimensional calculation having an additional momentum 
equation to solve plus the additional terras in all equations because of the 
presence of the third dimension. There is also additional overhead costs 
since the three-dimensional cases uses mass storage to transfer information 
in and out of core while the two-dimensional cases always have all the 
necessary information in core. The PEPSIS computer code was developed as a 
very general research tool and hence no attempt has been made to optimize its 
computational efficiency for a specific class of problems. The run times 
would decrease significantly if the number of terras in the governing 
equations could be reduced as for instance would occur for simplier geometric 
configurations. The present version has a general orthogonal capability 
as well as a limited nonor thogonal capability (see section on the 
User * s Manual ) . 
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V. DISCUSSION AND CONCLUSIONS 


The primary objective of this investigation was to develop and validate 
an efficient numerical procedure for the calculation of two- and 
three-dimensional supersonic flows (with embedded subsonic regions) in 
internal flow devices. It is felt that this objective has been achieved, and 
a new set of governing equations has been developed which are well-posed for 
solution by an efficient spatial forward marching procedure. This procedure 
has been validated by application to a series of test cases characteristic of 
the pehnomena that occur in internal flow devices and the results in general 
give very good agreement with the available experimental data. Numerous 
additional test cases have been successfully run by the present authors and 
by other investigators (Refs. 20 and 21) and branching or other unstable 
behavior has not been observed even when the marching step size has been 
several orders of magnitude below the minimum of the stability criterion that 
restrict schemes which have a minimum step size. Thus, the restrictive 
requirement of a minimum marching step size has not been observed for this 
set of governing equations. It should be emphasized that no approximation 
was made to the streamwise pressure gradient in the present approach. This 
is viewed as being extremely important for internal supersonic flows where so 
many of the phenomena of interest have associated with the large streamwise 
pressure gradients. Finally, the ef f iciency „ of the spatial forward marching 
procedure is such that a million grid point calculations can be performed in 
approximately four hours of IBM 370-3033 CPU run time (this is equivalent to 
approximately one hour of CRAY I run time). This procedure could be used on 
a routine basis for design type calculations for internal flow devices such 
as supersonic inlets. 
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IV. USER’S MANUAL 


The PEPSIS users' manual is meant to serve as a guide in helping the 
user make successful runs with the PEPSIS computer program. The degree of 
success obtained by the user will depend on the skill of the user and his 
ability to correctly apply the code to his particular problem. The code will 
solve the governing equations, subject to the user supplied boundary 
conditions, however, meaningful results will only be obtained if the boundary 
conditions are appropriate to the problem. In addition the user must specify 
viscosity models, initial conditions, a coordinate system and the location of 
grid points to adequately resolve the flow. The user with a good knowledge 
of the physics involved in his problem and how the code models the physics 
should, with a moderate amount of experience, be able to successfully apply 
the code to a wide variety of supersonic flow problems. 

The users' manual is divided into eight parts consisting of: (1) a flow 
diagram, (2) a brief description of each subroutine and its use, (3) a list 
of the Fortran variables and a description of their meaning, (4) a 
description of the logical file units utilized by the PEPSIS computer code, 

(5) a detailed description of the input required by the PEPSIS computer code, 

(6) a description of the common error conditions that may be encountered 
during the execution of a PEPSIS run and the corrective action to be taken, 

(7) sample input for two and three-dimensional cases and (8) sample output 
for the corresponding cases. 


Flow Diagram 

The purpose of the flow diagram is to help the user understand the basic 
flow of information within the PEPSIS computer code. Because of the size of 
the code (approximately 13,000 cards), a detailed flow diagram would be 
prohibitively large and probably be of little value to the user. Therefore, 
the flow diagram is intended only to give a general overview of the structure 
of the code. The interested user is urged to consult the program listing for 
details . 
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Flow Diagram for the PEPSIS Computer Code 
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PEPSIS Subrout ines 


Subroutine 

ADDRES 

ADI 

ADICP 

AD I UN 

AMARCH 

AMATRX 

ARTVIS 

AVRG 

BC 

BLKDATA 

BLT 

BULEEV 

CONVCT 

CORBND 

CORTRN 

CROSEC 

CURVT 

DATAS 

DELTX 

DELTXZ 


Purpose 

Calculate addresses for finite difference 
representation of metric and fluid dynamic 
variables . 

Master control subroutine for ADI procedure. 

Control subroutine for coupled equations. 

Control subroutine for uncoupled equations. 

Linearizes streamwise convective terms. 

Linearizes all streamwise terms. 

Artificial dissipation subroutine. 

Calculates averaged quantities in cross 
plane . 

Boundary condition subroutine. 

Stores default values of key variables. 

Calculates boundary layer thickness. 

Calculates Buleev turbulent mixing length. 

Linearizes cross plane convective terms. 

Calculates geometry transformation 
information on boundaries. 

Calculates geometry transformation 
information for interior points. 

Control subroutine for calculation of derived 
variables . 

Linearizes curvature terms. 

Logical file control subroutine. 

Calculates transformation information for 
ICORD = 2 option. 

Calculates transformation information for 
ICORD = 3 option. 
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Subroutine 

Purpose 

DIFF 

Linearizes diffusion terms* 

DISFCN 

Calculates dissipation function. 

DIV 

Calculates divergence of velocity. 

D0P2 

Control subroutine for linearization of 
Y-direction and source terms. 

DOP3 

Control subroutine for linearization of 
Z-direction terms • 

EOS 

Equation of state subroutine linearizes and 
updates pressure and temperature. 

FGFUN 

Calculates geometry groupings. 

GAUSS 

Solves uncoupled tri-diagonal set of 
equations. 

GENCBC 

Control subroutine for coupled boundary 
conditions. 

GENUBC 

Control subroutine for uncoupled boundary 
conditions • 

GEORD 

Controls reading of metric information from 
logical file unit LDRUM. 

GEOTRB 

Generates metric information on logical file 
unit LDRUM. 

INDIC 

Determines if flow is subsonic or supersonic 
at grid points. 

INPUTS 

Input subroutine. Input data enters and is 
processed. 

INTEBC 

Performs a two-dimensional linear 
interpolation for wall transpiration rates. 

LAMP 

Calculates laminar profile. 

LAW 

Calculates nondimens ional velocity, U + as a 
function of nondimens ional distance, Y + . 

LENGTH 

Calculates turbulent mixing length. 

LOADUP 

Brings into core dependent and derived 
variables at appropriate points. 
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Subroutine 


MAIN 

MATPRT 

MGAOSS 

MGERR 

NMLIST 

OUTPUT 

PLOT 

PLOTIN 

PROF 

QUICK 

READZ 

RESTRT 

ROTATE 

SETBVL 

SHEAR 

SPREAD 

SUB 

SWITCH 

TANHYP 

TNDER 


Purpose 

Main control program. 

Prints elements block tridiagonal matrix. 

Control subroutine for solving block 
tridiagonal systems of equations. 

Calculates error associated with solving 
block tridiagonal system of equations. 

Subroutine for printing namelist input 
Information. 

Control subroutine for printing out results 
on a cross-sectional (Y-Z) plane. 

Writes plot information on logical file unit 
JPLOT . 

Writes first record of general information on 
logical file unit JPLOT. 

Generates initial profiles. 

Matrix elimination subroutine. 

Prepares variables for printing. 

Reads and writes restart information. 

Rotates data from columns to rows and vice 
versa. 

Updates boundary information a line at a 
time . 

Control subroutine for the calculation of 
wall shear velocity. 

Spreads two-dimensional data to three 
dimensions . 

Contain special subsonic logic. 

Calculates streamwise location for switch of 
boundary condition. 

Grid stretch subroutine. 

Calculates normal derivative of temperature 
at a wall. 
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Subroutine 


Purpose 


TRANS 

TURB 

TURBP 

VISCOS 

WALLFN 

WRMATR 

YCALC 

ZERO 


Transition model subroutine. 

Turbulence model subroutine. 

Calculates turbulent profile based on theory 
of Maise-McDonald . 

Constant and laminar viscosity subroutine. 

Calculates wall shear velocity. 

Writes block tridiagonal dump information on 
logical file device NUNERR. 

Calculates Y and Z locations. 

Zeros out linearization arrays. 
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Logical File Units Utilized by PEPSIS Computer Code 


The PEPSIS COMPUTER code utilizes up to twelve (12) logical file units 
during the execution of a run stream. In many cases not all twelve units are 
used, and hence in these cases there is no need to define all twelve units. 
All references to a logical file unit in the PEPSIS computer code is 
accomplished through the use of a FORTRAN name rather than through a specific 
unit number. Thus, if the user desires to change a logical file unit number, 
this can be done through the input file. A list of the logical file units 
utilized by the PEPSIS computer code, their FORTRAN name, default value unit 
number, and a brief description of the use of the unit is presented below. 

All units are sequential. 


FORTRAN Name Default Unit Number 

MIN 5 

MOUT 6 

MASS1 8 

MASS2 9 


MSDD 15 

JDRUM 11 

LDRUM 12 

KDRUM 13 


Description 

Input data unit. 

Printed output unit. 

First unit which stores dependent 
and derived variables either by 
rows or columns. Not needed for 
two-dimensional cases, i.e., when 
TWOD = .TRUE. 

Second unit which stores dependent 
and derived variables either by 
rows or by columns . Not needed 
for two-dimensional cases, i.e., 
when TWOD = .TRUE. 

Unit which stores dependent and 
derived variables by rows only. 

Not needed for two-dimensional 
cases, i.e., when TWOD = .TRUE. 

Unit Which contains output of ADD 
computer code. Only needed when 
IGEOM = 10 or 11. 

Intermediate unit used in 
generations of final metric 
information file. Only needed 
when IGEOM = 10 or 1 1 . 

Unit which stores final metric 
information. Needed in all cases. 
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FORTRAN Name 


Default Unit Number 


Description 


NUNEER 


JPLOT 


JRSTIN 


JRSTOT 


14 Unit which stores information 

concerning the block tridiagonal 
matrix inversion. Needed when 
MGDMP * 0. 

16 Unit which stores plotting 

information. Needed when IPLOT 

* 0 . 

10 Input retart unit. This unit 

contains appropriate common block 
information and the value of the 
dependent and derived variables at 
each cross-sectional grid point at 
the restart streamwise station. 
Needed only when IRSTIN * 0. 

10 Output restart unit. This unit 

contains appropriate common block 
information, and the values of the 
dependent and derived variables at 
each cross-sectional grid point at 
the restart streamwise station. 
Needed only when IRSTOT = 0. 
Default is JRSTIN = JRSTOT; 
however, it is desired to have 
separate input restart and output 
restart files set JRSTOT = 17. 
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PEPS IS Input 


Except for an initial title card and plot file input data the entire 
PEPSIS input is entered by means of the NAMELIST format. There are two 
primary advantages to the use of the NAMELIST format: (1) if the default 

values (defined in the block data subroutine) are acceptable, the user need 
not input that variable, and (2) the order (within a given NAMELIST) in which 
the variables are entered is irrelevant. There are six NAMELIST input files 
in the PEPSIS code, $REST and $LIST1 through $LIST5. The first file is read 
in the main program and enters restart information. The remaining NAMELIST 
files are read in subroutine INPUTS. Basically, the NAMELISTS $LIST1 through 
$LIST5 can be divided by function. $LIST1 enters information about the 
governing equations and appropriate boundary conditions, $LIST2 enters 
reference and free stream conditions, $LIST3 enters geometric information, 
$LIST4 enters viscosity model and initial profile information and $LIST5 
enters file output information. A description of all the PEPSIS input 
information will be given below. 


Plot File Input 


Card 1 

Columns 

Format 

Variable 

Function 

1-24 

6A4 

TITLE(I) 

Title Card 

Card 2 

Columns 

Format 

Variable 

Function 


1-2 

112 

ISYM 

Reciprocal of Symmetry 

3-12 

1F10.0 

SYSTEM 

SYSTEM = 1 — Quasi- 
Cartesian Coordinates 




SYSTEM = 2 - Quasi- 

Cylindrical 

Coordinates 
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Namelist Input Description 


Namelist 

variable 

REST 

IRSTIN 

IRSTOT 

JRSTIN 

JRSTOT 

NFILE 
N SAVED 

ICOMP 


or 

name 


Description 
Restart Options 


Marching station number when data is to be read for 
restart case. 


IRSTIN = 0: Dead start case. 

IRSTIN t 0: Restart case started at station 

IRSTIN. 

Default value is 0. 

Interval for saving restart information. 

IRSTOT = 0: No restart information is saved. 

IRSTOT t 0: Information is saved at each IRSTOTth 

station. 

Default value is 0. 

Logical file name of input restart file. 

Default value is 10. 

Logical file name of ouput restart file. 

JRSTOT and JRSTIN do not have to be same file. 
Default value is 10. 

File number on unit JRSTIN desired for restart. 
Default value is 0. 

Number of restart stations saved on JRSTOT. 

On a restart by setting JRSTOT = JRSTIN and NFILE + 
NSAVED, one file can be used for both reading and 
writing without destroying the information 
previously saved. 

Default is -1. 


Flag for computer opt: 


ICOMP = 1 
ICOMP = 2 
ICOMP = 3 
ICOMP = 4 


Univac computer option. 

CDC computer option. 

IBM computer - virtual memory option. 
Disk writing computer option 


Default value is 4. 
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Namelist Input Description 


LIST1 

IHSTAG 


IBOUND(IBC) 


IEQ BC (IBC, 
JEQBC(IBC, 


Description 

Equations and Boundary Conditions 

IHSTAG * 0: Energy equation formulated in terms of 

static enthalpy, 

IHSTAG * 1: Energy equation formulated in terms of 

stagnation enthalpy, 

IHSTAG = 2: Stagnation enthalpy constant. 

Default value is 1. 

Computational domain boundary characteristics (wall 
or non-wall). 

IBOUND(IBC) = 1: Solid wall boundary at surface 

IBC (see list of FORTRAN 
variables for definition). 

IBOUND(IBC) * 2: Non-wall boundary at surface IBC. 

Default values are 1, 1, 2, 1. 

IEQ) Boundary condition of the governing equation IEQ at 

sold wall boundary IBC. 

Default values are: 12*2, 4*16, 4*11, 8*2. 

IEQ) Boundary condition of the governing equation IEQ at 

non-wall boundary IBC. 

Default values are: 4*11, 2*2, 4*11, 2*2, 16*11. 

Boundary condition options used either IEQBC or 
JEQBC are as follows: 

<f> : Any dependent variable 

P: Pressure 

T: Temperature 

n= Normal to boundary 
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Namelist Input Description 


Description LIST1 

Equations and Boundary Conditions 

subscript c* Cartesian component 


Index 

1 

2 

3 

4 

5 

6 

7 

8 


11 


12 

13 

14 


15 


Function Type 
Description 

A<}> * 0 (No change of $ at 
boundary) 

<t> = 0 

V or W known 
pV or pW known 
AP = 0 

P « PRESS(IBC) 

AT = 0 

T = TWALL(IBC) 

Derivative Type 
Description 

* 0 (gradient of $ normal 
o n 

to boundary) 

Mach line extrapolation using one- 
sided difference 

Slip boundary condition for velocity 
using wall function 

3 P 

= 0 (gradient of pressure normal 
to boundary) 

3P 

•tt— « curvature (pressure gradient 
o n 

normal to boundary with curvature 
effects) 



Namelist 

Input Description 



Derivative Type 

LIST1 

Index 

Description 


16 

Momentum equation in direction 
normal to boundary 


17 

= 0 (adiabatic condition for wall 
or symmetry condition for non-wall) 


18 

3T 

= DTDN(IBC) 


19 

Wall function boundary condition 
for temperature 


20 

3 A 

-j— * 0 (same as 11, but applied at 



one grid point off the wall) 


21 

Mach line extrapolation using central 
difference scheme at one point off the 
boundary 


41 

- 0 



3n 2 


42 

a2p = 0 



3n 2 


43 

o 



3n 2 


44 

3U 

3n c ° 


45 

3V 



0 


46 

n c *VP = 0 


47 

n c *VT = 0 
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LIST1 


Description 


TWALL(IBC) 


PRESS(IBC) 


DTDN(IBC) 


ASW(IBC) 

BSW(IBC) 

CSW(IBC) 

DSW(IBC) 


LIST2 

IUNITS 


LREF 
RE PL 
MINF 
PINF 


Specified temperature at boundary IBC. 

Default values are 4*1.0. 

Specified pressure at boundary IBC. 

Default values are 4*0.0. 

Specified temperature gradient at boundary IBC at 
boundary condition. 

Default values are 4*0.0. 

Coefficient of a cubic polynomial fit for IBCth 
surface to determine the axial location where 
boundary characteristics at boundary IBC should be 
changed from wall to non-wall or vice versa, i.e., 
IBOUND(IBC) automatically changed. 

Default values are: 

ASW = 4*1. 0E + 10 

BSW = 4*0.0 

CSW = 4*0.0 

DSW = 4*0.0 

Freest ream and Reference Conditions 


Sentinel for units. 

IUNITS * 1: English units 

IUNITS = 2: Metric units 

Default value is 2. 

Reference length (ft or m) 

No default. 

Reynolds number per unit length. 

No default. 

Free stream Mach number. 

No default. 

2 2 

Free stream static pressure (lbf/ft or nt/m ) 
No default . 
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LIST2 


Reference and Reference Conditions 


PZERO 


2 2 

Free stream stagnation pressure (lbf/ft or nt/m ) 
No default. 


PR 

PRT 

XOB(IBC) 


ICORD 


NE(IADI) 


NS 


Laminar Prandtl number. 

Default value if 0.74. 

Turbulent Prandtl number. 

Default value if 1.0. 

Location on IBCth surface where the boundary type 
switches from non-wall to wall. 

Default values are 4*1. 0E + 06. 

Flag for coordinate transformation. 

ICORD = 1: Conformal coordinates 

ICORD = 2: Nonorthogonal coordinates 

X + X 

Y + Y 

Z * C(X,Z) 

ICORD = 3: Nonorthogonal coordinates 

X + X 

Y ♦ n(X,Y,Z) 

Z + Z 

Default value is 1. 

Number of grid points in the Y(IADI - 1) and 
Z(IADI) = 2) directions. 

No default values. 

Number of last streamwise stations which solution 
is to marched to. 

No default value. 


61 


LIST2 


Freestream and Reference Conditions 


XENTR is the initial streamwise location. DELX is 
the initial stepsize in the streamwise (marching 
direction, i.e., X(2) - XENTR + DELX. At 
streamwise station I the streamwise position is 
given by X(I) = X(l-l) + AP(X(I-1) - X(I-2) where 
if AP is greater than 1.0, the streamwise step size 
will increase by (AP-1.0) percent each step. If AP 
is less than 1.0, the streamwise step size will 
decrease by (1.0-AP) percent each step DXMIN and 
DXMAX are lower and upper overriding limits on the 
step size. AP, DXMIN and DXMAX are dimensional so 
that streamwise step size variation can be changed 
by the IAP parameter, IAP denoting the streamwise 
location where these variables change. Values of 
XENTR and DELX should normally only be set on the 
initial run as these variables are automatically 
calculated for restarts. 

XENTR = 0.0 

DELX * No default value 

IAP = 1,9*1000000 

AP = 10*1.0 

DXMIN = 10*0.0 

DXMAX = 10*1. 0E + 06 

IGEOM Flag for coordinate options 

IGEOM = 1: Cartesian coordinates 

IGEOM = 2: Cylindrical coordinates 

IGEOM = 3: Polar coordinates 

IGEOM =10: General orthogonal coordinates 

(Cartesian in cross plane) 

IGEOM =11: General orthogonal coordinates 

(axi symmetric) 

Default value is 1 . 


XENTR, DELX, 

IAP(10), AP(10), 
DXMIN(IO), DXMAX(IO) 
(CONTINUED) 


XENTR, DELX, 

IAP(10), AP(10), 
DXMIN(IO), DXMAX(IO) 
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LIST3 


Geometric Options 


TWOD 


TT1 ( 2 ) 


TT2 (2 ) 


YS( 2 , 2) 


LIST4 

IBCP 


DELTAP(IBCP) 


Sentinel for two-dimensional option. 

If TWOD = .TRUE. TWO DIMENSIONAL 

TWOD « .FALSE. THREE DIMENSIONAL 

Default value is .FALSE. 

Grid distribution factor (lower surface (IBC = 1) 
Y-direction, left surface (IBC = 3) - 
Z-direc tion) . The closer the value is to 1.0, the 
tighter the packing; value must be negative. 

Default values are 2*0.0 

Grid distribution factor (upper surface (IBC = 2) 
Y-direction, right surface (IBC = 4) - 
Z-direct ion) . The closer the value is to the 1.0, 
the tighter the packing. 

Default values are 2*0.0 

Defines computational domain in Y-Z cross plane. 

YS(l,i) = 0.0 - lower limit Y-direction 

YS(2,1) “ 1.0 - upper limit Y-direction 

YS(1,2) = 0.0 - lower limit Z-direction 

YS(2,2) - 1.0 - upper limit Z-direction 

Default values are 0.0, 1.0, 0.0, 1.0. 

Initial Profile, Turbulence Information 


Basic surface for initial 

Boundary Layer Profile 

at surface 1 

at surface 2 

at surface 3 

at surface 4 
Default value is 1. 


profile generation. 


IBCP = 1 
IBCP = 2 
IBCP = 3 
IBCP = 4 


Boundary layer thickness on surface IBCP needed to 
generate the initial profile referenced to each 
surface . 


No default values. 
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LIST4 


Initial Profile, Turbulence Information 


CFP(IBCP) Skin friction coefficient on surface IBCP needed to 

generate the initial turbulent boundary layer 
profile . 

No default values, 

IPROF Flag for initial profile options. 

IPROF = 1: Freestream profiles 

IPROF = 2: Initial profiles supplied by user 

IPROF = 3: Boundary layer profiles based on 

necessary input 

IPROF = 4: Same as IPROF « 3, but angular 

components are obtained for general 
orthogonal coordinates 
Default value is 1. 

IMIXL Flag for mixing length options. 

IMIXL = 1: McDonald-Camarrata mixing length 

model based on prescribed boundary 
layer thickness (DELTAB) . With wall 
shear value used to calculate non- 
dimensional distance 

IMIXL = 2: Buleev mixing length model 

IMIXL = 3: McDonald-Camarrata mixing length 

model based on dynamically obtained 
boundary layer thickness with fixed 
wall shear 

IMIXL = 4: Same as IMIXL = 1, but local shear is 

used to calculate nondimens ional 
distance 

IMIXL = 5: Same as IMIXL = 3, but local shear 

is used 
Default value is 1 . 

BETA Angle of attack in degrees. 

Default value is 0.0. 

YAW Yaw angle in degrees. 

Default value is 0.0. 


64 


LIST4 


Initial Profile, Turbulence Information 


IVISC Flag for viscosity options. 

IVISC = 1: Constant viscosity 

IVISC = 2: Laminar viscosity obtained from 

Sutherland 's relation 

IVISC = 3: Turbulent viscosity is obtained 

from mixing length model 
Sutherland's law for laminar viscosity; 

IVISC = 4: Turbulent viscosity obtained from 

TKE - mixing length model 
Sutherland's law for laminar viscosity; 
Default value is 1. 

DELTAB(IBC) Specified boundary layer thickness on surface IBC 

for mixing length model of turbulence. 

No default values. 

ITRANS Flag which tells whether transition turbulence 

model logic is used. 

ITRANS = 0: No transitional model is used 

ITRANS * 0: Transitional model is used 

Default value is 0. 

IBLT Flag which tells whether boundary layer thickness 

is input or calculated dynamically. 

IBLT = 0: Boundary layer thickness is input 

IBLT * 0: Boundary layer thickness is dynamically 
calculated. 

Default values is 0. 

TKEINF Freestream turbulent kinetic energy. 

is 0.0. 

Default value is 0.0. 

IVARPR(I) Index of variables to be printed. Needed only for 

three-dimens ional flow. 

IVARPR(I) =0: No print 

IVARPR(I) = 1: Print every IPRINT steps 
IVARPR(I) = 2: Print every JPRINT steps 
Default values are 5*1, 2*0, 3+1, 7+0. 
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LIST4 


Initial Profile, Turbulence Information 


1=1: UVEL 

1=2: VUEL 

1=3: WVEL 

1=4: Density 

1=5: Enthalpy 

1=6: Turbulent kinetic energy 

1=7: Turbulent dissipation 

1=8: Pressure 

1=9: Temperature 

I = 10: Mach Nunber 
I = 11: Mach Number Indicator 
I = 12: Stagnation temperature 
I = 13: Stagnation pressure 
I = 14: Pressure coefficient 
I = 15: Laminar viscosity 
I = 16: Mixing length 
I = 17: Turbulent viscosity 
I = 18: Effective viscosity 

IPLOT Marching station interval for storage of plotting 

information • 

IPLOT =0: No plotting 

IPLOT * 0: Store plotting information every 

IPLOT station; IPLOT cannot be changed 
during a run. 

Default value is 0. 

IPRINT Primary marching station interval for printing. 

Default value is 1. 

JPRINT Secondary marching station interval for printing. 

Default value is 1. 


66 



Error Conditions in the PEPSIS Computer Code 


Failure of the PEPSIS computer code to successfully execute a runstream 
can occur because of either inconsistent or incorrect input data or because 
of an attempt to apply the PEPSIS code to a case where the physics violate 
the assumptions inherent in the code. This section will address only the 
former mode of failure. Avoidance of the latter failure mode is dependent 
primarily on the users understanding of the basic physics of the case he is 
going to run, and the degree to which the PEPSIS code can be expected to 
model the physics. 

One method of discussing the inconsistent or incorrect input data mode 
of failure is by examining the possible failures in the various subroutines. 
Since the individual subroutines are responsible for separate tasks during 
the execution of a run, (e.g. overall control of the program geometry 
generation, etc.), this technique will in essence outline the possible 
failure modes as the tasks are performed. Discussion will occur in the same 
order as the run is executed. 

SUBROUTINE RESTRT 

There are two modes by which SUBROUTINE RESTRT can fail. Both involve 
improper use of the restart file. A message, PESTART INFORMATION REQUESTED 
AT (IRSTIN marching number) BUT STORED INFORMATION AT SEQUENCE (NFILE) IS AT 
STATION (Station number). This message occurs because the marching station 
number read off the NFILEth restart file does not match the input value of 
IRSTIN. The corrective action is to make NFILE and IRSTIN consistent with 
each other. Another possible mode of failure occurs when NFILE exceeds the 
number of files on the restart device, JRSTIN, in which case an END of 
INFORMATION (or analogous statement) will appear in the day file. The 
corrective action is to recheck the input value of NFILE. If JRSTIN * 

JRSTOT, the value of NFILE is the number of the restart on device JRSTIN. 

SUBROUTINE INPUTS 

There are two failure modes in SUBROUTINE INPUTS. In the first case, 
the message NS = (input value of NS) GREATER THAN NSMAX = (dimension of X 
vector) will be printed if the number of marching stations exceeds the 
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dimensioned value of X, the streamwise locations. The corrective action is 
to lower the value of NS. The second failure mode occurs when the Buleev 
turbulence model is specified for a two-dimensional case. Since this model 
is not applicable to two-dimensional cases, the message CANNOT USE BULEEV 
TURBULENCE MODEL IN TWO-DIMENSIONAL FLOW is printed. The corrective action 
is to specify an alternate turbulence model. 

SUBROUTINE GEOTRB 

At present SUBROUTINE GEOTRB is coded to calculate metric information 
for values of IGEOM = 1, 2, 3, 10 and 11. Values of IGEOM 4-9 are left for 
various coordinates that may be coded in the future. Input value of IGEOM = 
4-9 will result in the message INVALID OPTION IN GEOTRB. The corrective 
action is to either change the value of IGEOM or to code in a new option. 

For IGEOM options 10 and 11 (conformal-Cartesian cross-section and 
conformal-axisymmetric cross-section) the metric information is externally 
generated by the ADD computer code. In this case, logical file units JDRUM 
and KDRUM must be defined. JDRUM contains the ADD code data which is then 
interpolated onto the PEPSIS mesh system. If the PEPSIS values of the 
streamwise coordinate is less than the first value of the ADD code streamwise 
coordinate no streamwise interpolation is possible and the message FAILURE IN 
GEOTRB - SQ12 = (PEPSIS position) SQ1 = (first ADD code position) SQ2 = 
(second ADD code position). The corrective action is to increase the value 
of XENTR (the first PEPSIS position) to a value greater than SQ1 . On the 
other hand, if the value of a PEPSIS streamwise coordinate exceeds the last 
streamwise position generated by the ADD code an END OF INFORMATION message 
will appear in the day file. The corrective action is to either rerun the 
ADD code such that the maximum PEPSIS streamwise coordinate does not exceed 
the maximum ADD code streamwise coordinate or to reduce the maximum PEPSIS 
streamwise coordinate to an acceptable value. 

SUBROUTINE INTEBC 

SUBROUTINE INTEBC performs a two-dimensional linear interpolation of 
the transpiration schedules on both X-Y planes at X-Z planes. If the number 
of streamwise stations on a surface at which data is input exceeds 15 (the 
dimensioned size of the data arrays) a message FAILURE IN INTEBC VALUE OF 
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NPTSX (surface number IBC) « (value of NPTSX(IBC) EXCEEDS DIMENSION LIMITS OF 
15 is printed. The corrective action is either to updimension NPTSX and 
associated variables or to decrease the value of NPTSX. Likewise, in the Y 
or Z direction data can be input at up to 15 locations. If the value of 
NPTSYZ exceeds 15, the message FAILURE IN INTEBC VALUE OF NPTSYZ {streamwise 
location, surface number) = (value of NPTSYZ) EXCEEDS DIMENSION LIMITS OF 15 
Is printed. The corrective action is either to updimension NPTSYZ and 
associated variables or to decrease the value of NPTSYZ. 

SUBROUTINE QUICK 

If the choice of boundary conditions is incorrectly made, it is 
possible that a singular matrix will result. This will manifest itself in 
SUBROUTINE QUICK in an attempt to divide by zero. The corrective action is 
to re-evaluate the choice of input boundary conditions to determine the 
source of the singularity. An example of an improper choice of a boundary 
condition set would be to choose as boundary conditions the three no-slip 
conditions for the three momenta equations, the normal pressure condition for 
the continuity equation and the normal momentum equation for the enthalpy 
equation. In this case, the enthalpy does not appear in any of the boundary 
conditions, and hence a singular matrix would result. 

SUBROUTINE CROSEC 

Often, if a case is not going to successfully run, the code will cease 
operation in SUBROUTINE CROSEC. This will occur because of the existence of 
a negative temperature in which case the Mach number calculation will fail in 
SQRT. There can be many reasons for this failure mode. Usually, however, it 
can be related to inadequate numerical resolution of the physical processes 
that are occurring. For instance, a lack of transverse grid points might 
lead to large oscillations in the pressure or too large a streamwise step in 
the region where a wall inclination is rapidly changing might result in a 
temperature becoming negative. Sometimes it is difficult to know a priori 
what grid resolution is necessary for a given problem. Usually, 
experimentation with two-dimensional cases can provide some guidelines for 
three-dimensional cases. This in addition with the users' overall experience 
with the code and his understanding of the physical processes will usually 
provide the means of resolving the above problem. 
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PEPSIS FORTRAN VARIABLES 

FORTRAN 

COMMON 

DESCRIPTION 

SYMBOL 

BLOCK 


ACON 

LAUU 

CONSTANT IN ARGUMENT OF EXPONENTIAL FUNCTION 
FOR TRANSITIONAL MODEL 

AG(NN,9,2) 

OPER 

DIFFERENCE WEIGHTS IN PHYSICAL COORDINATES 

AGEO 

GEOM 

COEFFICIENTS OF POLYNOMIAL FIT FOR BOUNDARY SHAPE 

AG1D(9) 

FGCOM 

TEMPORARY STORAGE ARRAY OF METRIC INFORMATION 

AG1P 

FGCOM 

TEMPORARY STORAGE ARRAY OF METRIC INFORMATION 

AG2D( 9) 

FGCOM 

TEMPORARY STORAGE ARRAY OF METRIC INFORMATION 

AHP(5) 

FGCOM 

TEMPORARY STORAGE ARRAY OF METRIC INFORMATION 

AH1D<5,9) 

FGCOM 

TEMPORARY STORAGE ARRAY OF METRIC INFORMATION 

AIE(NN f 7) 

PRFILE 

INITIAL PROFILE ARRAY 

AM(NCPLD, 3ANCPLD+1 )CCOH 

UTILITY MATRIX USED IN BLOCK MATRIX INVERSION 

AMACRT 

SUPER 

MACH NUMBER CRITERION USED IN LOCATING SONIC LINE 

AN<NEQS,NN) 

LIN 

STORAGE FOR LINEARIZATION COEFFICIENTS OF 
X - DERIVATIVES 

AP(10) 

GEOM 

AMPLIFICATION RATE OF MARCHING STEP SIZE 

APLUS 

LAUU 

CONSTANT IN ARGUMENI OF EXPONENTIAL 
FUNCTION FOR VAN DRIEST DAMPING FORMULA 

ASU(4) 

BOUND 

COEFFICIENTS OF POLYNOMIAL FIT FOR SWITCHING THE 
BOUNDARY SURFACE IYPE 

AVISC(2,NEQS) 

vise 

COEFFICIENT USED IN ARTIFICIAL DAMPING 

BETA 

REF 

ANGLE OF ATTACK 

BGEO 

GEOM 

COEFFICIENTS OF POLYNOMIAL FIT FOR BOUNDARY SHAPE 

BL0CK1 ( IADD3) 


STORAGE FOR ADD CODE MEIRIC INFORMATION 
EQUIVALENCED TO C( 1,1,1) 

BL0CK2( IABD3) 


STORAGE FOR ADD CODE METRIC INFORMATION 
EQUIVALENCED TO C< 1 , 1 ,NN/2+l ) 
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FORTRAN 

COMMON 

DESCRIPTION 

SYMBOL 

BLOCK 


BLTH(NN f 4) 

LAUU 

DYNAMICALLY DETERMINED BOUNDARY LAYER THICKNESS 

BSU(4) 

BOUND 

COEFFICIENTS OF POLYNOMIAL FIT FOR SNITCHING THE 
BOUNDARY SURFACE TYPE 

BUD 

LIN 

CRANK-NICHOLSON FACTOR 

BUD I 

LIN 

INVERSE OF BUD 

C(NN,NCPLD,NN) 

CCOM 

BLOCK DIAGONAL MAIRIX ELEMENIS 

CDUM(NDIM) 


TEMPORARY SIORAGE ARRAY TO ROTATE DATA FROM COLUMNS 
TO ROUS AND VICE VERSA - EQUIVALENCED TO C( 1,1,1) 
NDIM = MZVAR A MLEVEL A NN 

CFP(4) 

PRFILE 

SKIN FRICTION COEFFICIENT 

CGEO 

GEOH 

COEFFICIENTS OF POLYNOMIAL FIT FOR BOUNDARY SHAPE 

CMUINF 

vise 

C0N8TANT IN TURBULENT VISCOSITY MODEL 

CONGEO( 1 1 ANN ) 

GEOM 

COORDINATE TRANSFORMATION INFORMATION 

CONVDR 

UNITS 

CONVERSION FACTOR IN GOING FROM DEGREES TO RADIANS 

CONVRD 

UNITS 

CONVERSION FACTOR IN GOING FROM RADIANS 10 DEGREES 

COOR(NN,4) 

GEOM 

PHYSICAL COORDINATE INFORMATION UITH RESPECT 
TO ABSOLUTE ORIGIN AT N+1ST STREAMUISE LOCATION 

C00RN(NN,4) 

GEOM 

PHYSICAL COORDINATE INFORMATION UITH RESPECT 
TO ABSOLUTE ORGIN AT NTH STREAMUISE STATION 

CPINF 

FREE 

FREE STREAM SPECIFIC HEAT 

CPREF 

REF 

REFERENCE PRESSURE COEFFICIENT 

CPREFI 

REF 

INVERSE OF REFERENCE PRESSURE COEFFICIENT 

CRITU 

OPER 

CRITICAL VELOCITY USED FOR FLARE APPROXIMATION 

CSOLNtNCPLD , NN) 

CCOM 

SOLUTION 10 BLOCK TRI-DIAGONAL MATRIX INVERSION 

CSU ( 4) 

BOUND 

COEFFICIENTS OF POLYNOMIAL FIT FOR SUIICHING THE 
BOUNDARY SURFACE TYPE 
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FORTRAN 

SYMBOL 

COMMON 

BLOCK 

DESCRIPTION 

CTWO 

REF 

NUMERICAL CONSTANT IN GOVERNING EQUATION 

CXI(NDIFh) 

ADDR 

STORAGE FOR FIRST DERIVATIVE DIFFERENCE WEIGHTS 

CXXI(NBIFM) 

ADDR 

STORAGE FOR SECOND DERIVATIVE DIFFERENCE WEIGHTS 

Cl(NN) 

CCOM 

SUBDIAGONAL MATRIX ELEMENTS 

C1SUIH 

vise 

COEFFICIENT IN SUTHERLAND'S LAW OF LAMINAR 
VISCOSITY 

C2(NN) 

CCOM 

DIAGONAL MATRIX ELEMENTS 

C2SUTH 

vise 

COEFFICIENT IN SUTHERLAND'S LAW OF LAMINAR 
VISCOSITY 

C3(NN) 

CCOM 

SUPERDIAGONAL MATRIX ELEMENTS 

C4(NN) 

CCOM 

VECTOR ELEMENTS 

D 

VAR 

INDEX FOR DIVERGENCE 

DELMAX(NEQS) 

EQN 

MAXIMIUM VALUES OE THE DELTAS 

DELTAB< 4) 

LAUD 

SPECIFIED BOUNDARY LAYER THICKNESS FOR 
MIXING LENGTH MODEL OF TURBULENCE 

DELTAP (4) 

PRFILE 

BOUNDARY LAYER THICKNESS 

BELX 

GEOM 

STEP SIZE IN MARCHING DIRECTION 

BGEO 

GEOM 

COEFFICIENTS OF POLYNOMIAL FIT FOR BOUNDARY SHAPE 

BIFOP(6) 

OPER 

DIFFERENCE WEIGHTS IN COMPUTATIONAL COORDINATES 

Dt 

VAR 

STORAGE LEVEL OF DIVERGENCE IN 
GENERAL PURPOSE STORAGE 

BM1 (NCPLD, NCPLD) 

CCOM 

GENERAL PURPOSE STORAGE ARRAYS FOR BLOCK 
TRI-DIAGONAL MATRIX INVERSION 

DM2(NCPLD) 

CCOM 

GENERAL PURPOSE STORAGE ARRAYS FOR BLOCK 
TRI-DIAGONAL MATRIX INVERSION 

BM3 ( NCPLD , NCPLD ) 

CCOM 

GENERAL PURPOSE STORAGE ARRAYS FOR BLOCK 
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FORIRAN COMMON DESCRIPTION 

SYMBOL BLOCK 

TRI-DIAGONAL MATRIX INVERSION 

DS VAR INDEX FOR DISSIPATION FUNCTION 

DSL VAR STORAGE LEVEL OF DISSIPATION FUNCTION IN 

GENERAL PURPOSE STORAGE 

DSW(4) BOUND COEFFICIENTS OF POLYNOMIAL FIT FOR SWITCHING THE 

BOUNDARY SURFACE TYPE 

DTDNI4) BOUND SPECIFIED TEMPERATURE GRADIENT NORMAL TO 

BOUNDARY 


DIDNU(NN,4) 

BOUND 

STORAGE ARRAY FOR TEMPERATURE GRADIENT 
TO BOUNDARY 

NORMAL 

DX 

OPER 

STEP SIZE IN X-DIRECT ION 


DXI 

OPER 

INVERSE OF DX 


DXMAX(IO) 

GEOM 

MAXIMUM MARCHING STEP SIZE 


DXMIN(IO) 

GEOM 

MINIMUM MARCHING STEP SIZE 


D1(NEQS,NDIFM,NN) 

LIN 

STORAGE FOR LINEARIZATION COEFFICIENTS 
Y - DERIVATIVES 

OF 

D1L(NEQS,NDIFM) 

LIN 

STORAGE OF LINEARIZATION COEFFICIENT 
FOR TEMPERATURE AND PRESSURE COMPUTATION 

D2(NEQS,NDIFM r NN) 

LIN 

STORAGE FOR LINEARIZATION COEFFICIENTS 

OF 


Z - DERIVATIVES 

E(NN,NCPLD,NN) ARRAY USED IN MGAUSS ERROR CHECK - EQUIVALENCED 

TO C(l f l,l) 

EPS VAR INDEX FOR DISSIPATION OF TURBULENCE 

KINETIC ENERGY 

EPSMUF LAUW CONVERGENCE CRITERION FOR WALL FUNCTION FORMULATION 

F(14 r NN) ZPLOT TEMPORARY STORAGE FOR PLOT INFORMATION 
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FORTRAN 

SYMBOL 

COMMON 

BLOCK 

DESCRIPTION 

FACLM<4) 

LAUU 

MULTIPLICATION FACTOR TO BOUNDARY LAYER 



THICKNESS 

FG(8,3,NN) 

METRIC 

STORAGE FOR METRIC COEFFICIENTS AND DERIVATIVES 

GAM(3) 

GEOM 

COORDINATE TRANSFORMATION COEFFICIENT FOR 



BOUNDARY CONDITIONS 

GAMMA 

REF 

RATIO OF SPECIFIC HEATS 

GC 

UNITS 

GRAVITY CONSTANT 

H 

VAR 

INDEX FOR ENTHALPY 

HFORM 

REF 

HEAT OF FORMATION 

HINF 

FREE 

FREE STREAM ENTHALPY 

HREF 

REF 

REFERENCE ENTHALPY 

HREFI 

REF 

INVERSE OF REFERENCE ENTHALPY 

IA 

CCOM 

INDEX REFERRING TO SUBDIAGONAL MAIRIX ELEMENIS 

IADDO(NBIFM) 

ADDR 

ADDRESSES IN THE OPPOSITE DIRECTION 

IADDP (NBIFM ) 

ABDR 

ADDRESSES IN THE PRIMARY DIRECTION 

IADBSl(NN) 

ABBRF 

ADDRESS FOR POINT LOGIC OF FLUID VARIABLES 

IADDS2<NDIFM, NN) 

ADDRF 

ADDRESS FOR Y DERIVATIVE OF FLUID VARIABLES 

IADDS3(NDIFM,NN) 

ADBRF 

ADDRESS FOR Z DERIVATIVE OF FLUID VARIABLES 

IADDS4(NBIFMAA2, 

NN) ADDRF 

ADDRESS FOR MIXED DERIVATIVE OF FLUID VARIABLES 

IADDS5(NDIFMAA2, 

NN) ADDRF 

ADDRESS FOR MIXED DERIVATIVE OF FLUID VARIABLES 

IADD1 

ADD 

NO. OF GEOMETRIC VARIABLES USED IN ADD CODE 

IADD2 

ADD 

NO. OF TRANSVERSE GRIDPOINT USED IN ADD CODE 

IADD3 

ADD 

RECORD SIZE USED IN ADD CODE 

IADI 

SUEEP 

ADI SWEEP DIRECTION 

IADIM1 

SWEEP 

IADI - 1 

IAP ( 10) 

GEOM 

MARCHING STEP INDEX AT WHICH AP , DXM IN , DXMAX 
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FORTRAN 

SYMBOL 

IB 

IBC 

IBCP 

IBLT 

IBOUND (4) 

IC 

ICDC(NN+1 f 2) 
ICOMP 

i 

' ICONS (3, NEQS ) 

ICORD 

ICPLD(NEQS,2) 

ID 

ID IF( 40 ) 

IDMPY 

IDHPZ 

IDUM2(NN) 

IDUH3(NN) 

IEQ 

, IEQBC(4,NEQS) 


COMMON DESCRIPTION 

BLOCK 

ARE REINITIALIZED 

CCOM INDEX REFERRING TO DIAGONAL MATRIX ELEMENIS 

BOUND INDEX FOR BOUNDARY SURFACE IDENTIFICATION 

PRFILE BASIC SURFACE FOR INIIIAL PROFILE GENERATION 
LAWW FLAG UHICH TELLS WHETHER BOUNDARY LAYER 

THICKNESS IS INPUT OR CALCULATED DYNAMICALLY 
BOUND BOUNDARY SURFACE TYPE INDICATOR 

CCOM INDEX REFERRING TO SUPERDIAGONAL MAIRIX ELEMENTS 

CDC RECORD INDEX FOR READMS AND WRITEMS MASS 

STORAGE DEV ICES-CDC COMPUTER ONLY 
REST FLAG FOR COMPUTER OPTIONS 

OPER FLAG FOR CONVECIIVE FORMULATION BASED ON MACH 

NUMBER AND EQUATION 

GEOM FLAG FOR COORDINATE TRANSFORMATION OPTIONS 

EQN COUPLED EQUATION SENTINEL 

CCOM INDEX REFERRING TO VECTOR MATRIX ELEMENTS 

OPER INDEX FOR TYPE OF DIFFERENCING OF BOUNDARY 

CONDITIONS 

DMP DUMP LINE NO. IN Y DIRECTION 

DMP DUMP LINE NO. IN Z DIRECTION 

ADDRG INDICATOR OF THE POINT IN DIFFERENCE MOLECULE 

WHERE Y DERIVATIVE IS TAKEN 

ADDRG INDICATOR OF IHE POINT IN DIFFERENCE MOLECULE 
WHERE Z DERIVATIVE IS TAKEN 
EQN EQUATION NUMBER INDEX 

BOUND SPECIFIED BOUNDARY CONDITION FOR WALL 
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FORTRAN 

SYMBOL 

COMMON 

BLOCK 

DESCRIPTION 

IEQNUM ( NEQS ) 

EQN 

IDENTIFICATION NUMBER OF EQUATIONS TO BE SOLVED 

IFLARE 

OPER 

SENTINEL WHICH TELLS IF FLARE OPTION IS USED 

IGDMP 

BMP 

FLAG FOR DUMP OPTIONS 

IGEOM 

GEOM 

FLAG FOR COORDINATE SYSTEM OPIIONS 

IH 

GRID 

UPPER BOUNDARY POINT ON THE LINE WHERE IMPLICIT 
SOLUIION IS OBTAINED 

IHSTAG 

EQN 

FLAG FOR ENERGY EQUATION OPTIONS 

IJ 

DIFCOM 

COLUMN OR ROW NUMBER ON WHICH CALCULATION IS BEING 
MADE 

IL 

GRID 

LOWER BOUNDARY POINT ON THE LINE WHERE IMPLICIT 
SOLUTION IS OBTAINED 

IMIXL 

LAUU 

ELAG FOR MIXING LENGTH OPTIONS 

IND 

VAR 

FLAG WHICH TELLS IF MESH POINT BELONGS TO 
SUPERSONIC OR SUBSONIC REGION 

INDC(NN) 

DIFCOM 

MACH NUMBER INDICATOR 

INDL 

VAR 

STORAGE OF MACH NUMBER INDICATOR IN 
GENERAL PURPOSE STORAGE 

INH12 

ADD 

SENTINEL USED TO DETERMINE FORMULATION USED 
IN CALCULATION OF TRANSVERSE DERIVATIVE OF HI 

INH21 

ADD 

SENTINEL USED TO DETERMINE FORMULATION USED 
IN CALCULATION OF STREAMUISE DERIVATIVE OF H2 

INH31 

ADD 

SENTINEL USED TO DETERMINE FORMULATION USED 
IN CALCULATION OF STREAMWISE DERIVATIVE OF H3 

INH32 

ADD 

SENTINEL USED TO DETERMINE FORMULATION USED 
IN CALCULATION OF TRANSVERSE DERIVATIVE OF H3 

IOPTUF 

LAUU 

SENTINEL WHICH DETERMINES WALL EUNCTION FORMULATION 
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FORTRAN COMMON 

SYMBOL BLOCK 

IOPIYZ(3,NEQS,2) OPER 


IPLOI 

IPRINT 

IPROF 

IPRTE 

IRSTIN 

IRSTOT 


CIO 

CIO 

PRFILE 

GASL 

REST 

REST 


ISONIC SUPER 

ISS(NN,4) SUPER 

ISSHET SUPER 

ISTART HETRIC 

ISU UNIVAC 

ISYM ZPLOT 

ITRANS LAUU 

IUNITS UN ITS 

IVARNO(NEQS) EQN 

IVARPR ( 25 ) PRNT 

IVISC vise 

IUALE VISC 


DESCRIPTION 

FLAG FOR DIFFERENCING FORMULATION BASED ON 
MACH NUMBER, EQUATION, AND ADI DIRECTION 
MARCHING STATION INTERVAL FOR STORAGE OF 
PLOTTING INFORMATION 

PRIMARY MARCHING STAIION INTERVAL FOR PRINTING 

FLAG FOR INITIAL PROFILE OPTIONS 

FLAG WHICH DETERMINES EQ. OF STAIE FORMULATION 

STREAMWISE STATION NUMBER FOR RESTART 

STREAMUISE INTERVAL FOR SAVING RESTART 

INFORMATION 

FLAG FOR SONIC LINE INTERPOLATION LOGIC 

GRID POINT LOCATION OF SONIC LINE 

INDEX USED 10 DETERMINE IF SONIC LINE IS LAST 

SUBSONIC POINT OR FIRSI SUPERSONIC POINT 

INITIAL CONDITION INDEX 

SENTINEL FOR WORD ADDRESSABLE OR SECTOR- 

ORIENTED MASS STORAGE DEVICE - UNIVAC ONLY 

SYMMETRY OPTION FOR PLOTS 

FLAG WHICH TELLS WHETHER TRANSITION TURBULENCE 
MODEL LOGIC IS USED 

FLAG USED TO DETERMINE SET OF DIMENSION UNITS USED 

IDENTIFICATION NUMBER OF DEPENDENT VARIABLES 

INDEX OF VARIABLES TO BE PRINTED 

FLAG FOR VISCOSITY OPTIONS 

SENTINEL WHICH DETERMINES IF WALL FUNCTION 

LOGIC IS NEEDED IN THE CALCULATION OF WALL 
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FORTRAN 

SYMBOL 

COMMON 

BLOCK 

DESCRIPTION 



VISCOSITY 

IWR 

CIO 

SENTINEL FOR NAMELIST REST PRINT 

IYGD(NDIFM) 

ADDR 

GEOMETRY ADDRESSES 

IZCK3) 

ZEX1 

RELATIVE UNIT NO. FOR VIRTUAL MEMORY STORAGE 

11 IG(8 f 3) 

FGCOM 

INDEX NEEDED IN THE CALCULATION OF METRIC 



INFORMATION 

I2IG(8,3) 

FGCOh 

INDEX NEEDE IN THE CALCULATION OF METRIC 



INFORMATION 

J A < 3 ) 

ADDRF 

SHIFT LOGIC INDEX 

JADDO(NDIFM) 

ADDR 

ADDRESSES IN THE OPPOSITE DIRECTION 

JADDP(NBIFM) 

ADDR 

ADDRESSES IN THE PRIMARY DIRECTION 

JBOUNB(NN r 4) 

BOUND 

BOUNDARY TYPE INDICATOR AI EACH POINT ON BOUNDARY 

JDMAX 

EQN 

Y GRID POINT LOCATION OF MAXIMIUM DELIA 

JORUM 

CIO, 

ADD CODE DEVICE 

JDUM 

DIFCOM 

INDEX DENOTING RELATIVE POINT ABOUT WHICH 



DERIVATIVE IS LOCATED 

JEQBC(4,NEQS) 

BOUND 

SPECIFIED BOUNDARY CONDITION FOR NON-UALL 

JEQN(NEQS f 2) 

EQN 

EXTERNAL EQUATION NUMBER 

JGSTOR 

LIN 

VALUE OF JG NEEDED BY SUBROUTINE EOS 

JPLOT 

CIO 

DEVICE FOR PLOTTING 

JPR INI 

CIO 

SECONDARY MARCHING STATION INTERVAL FOR PRINTING 

JPROF ( 4 ) 

PRFILE 

SENTINEL FOR BOUNDARY VALUES DURING INITIAL PROFILE 



GENERATION 

JRSIIN 

REST 

LOGICAL FILE FROM WHICH RESTART INFORMATION 



IS READ 

JRSTOI 

REST 

LOGICAL FILE ON WHICH RESTART INFORMATION 
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FORTRAN COMMON DESCRIPTION 

SYMBOL BLOCK 

IS WRITTEN 

JVAR(NEQS,2) EQN VARIABLE NUMBER ASSOCIATED WITH EACH EQUATION 

DURING AN ADI SWEEP 

CIO SENTINEL FOR NAMELIST PRINT OPIION 

OPER RELATIVE MARCHING STATION COUNIER 

OPER ABSOLUTE MARCHING STATION COUNTER 

BMP MARCHING STAIION WHEN DUMP OUTPUT IS REQUESTED 

ADDRG INDICIES NECESSARY TO CALCULATE GEOMETRIC GROUPINGS 

EQN Z GRID POINT LOCAIION OF MAXIMIUM DELIA 

CIO DEVICE FOR TEMPORARY STORAGE - USED IN GEOMETRY 

GENERATION 

BOUND ADDRESSES FOR BOUNDARY CONDITIONS 

CIO DEVICE FOR FINAL METRIC INFORMATION 

EQN LOUEST INDEX OF EQUATIONS SOLVED EITHER BY 

COUPLED OR UNCOUPLED ADI SWEEP 

' LEQ2 EQN HIGHEST INDEX OF EQUATIONS SOLVED EITHER BY 

I 

COUPLED OR UNCOUPLED ADI SWEEP 
ADDR GEOMETRY LEVEL 

ADDRG GEOMETRY LEVEL 

ADDRG ADDRESS FOR POINT LOGIC OF GEOMETRIC VARIABLES 

ADDRG ADDRESS EOR Y - DERIVATIVE OF GEOMETRIC 
VARIABLES 

LGA3(NBIFM,NN) ADDRG ADDRESS FOR Z - DERIVATIVE LOGIC OF GEOMETRIC 

VARIABLES 

LGA4(NDIFMAA2,NN> ADDRG ADDRESS EOR CROSS DERIVATIVE( Y-Z) LOGIC OF 

GEOMETRIC VARIABLES 


LEV(3) 

LEVEL 

LGAKNN) 

LGA2(NDIFM.NN) 


LADD (3) 

LDRUM 

LEQI 


JWR(5) 

JX 

JXDUM 

JXDUMP 

KA<5) 

KDMAX 

KDRUM 
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FORTRAN 

SYMBOL 

LGA5<NDIFMAA2 

LREF 

LREFI 

LSHFT 

LVG(8) 

MASS1 

MASS2 

MCPLD 

MEFE 

MEFFL 

MEQK 

MEQS 

MEQS1 

MEQS2 

MGDMP 

MGD1 

MGB2 

MIN 

MINE 

ML 

MLEVEL 

MLL 

MN 
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COMMON 

BLOCK 

NN) ADBRG 

REF 

REF 

GEOM 

BOUND 

CIO 

CIO 

EON 

VAR 

VAR 

EQN 

EQN 

EQN 

EQN 

DMP 

GRID 

GRID 

CIO 

FREE 

VAR 

PARAM 

VAR 

VAR 


DESCRIPTION 

ADDRESS FOR CROSS DERIVATIVE(Z-Y) LOGIC OF 

GEOMETRIC VARIABLES 

REFERENCE LENGTH 

INVERSE OF REFERENCE LENGTH 

SHIFT INDEX FOR COORDINATE TRANSFORMATION 

BOUNDARY POINT INDICATOR FOR BOUNDARY CONDITION 

GENERAL PURPOSE HASS STORAGE DEVICE 

GENERAL PURPOSE MASS STORAGE DEVICE 

NUMBER OF COUPLED EQUATIONS TO BE SOLVED 

INDEX FOR EFFECTIVE VISCOSITY 

STORAGE LEVEL OF EFFECTIVE VISCOSITY IN 

GENERAL PURPOSE STORAGE 

LEQ1 - 1 

TOTAL NUMBER OF EQUATIONS TO BE SOLVED 
INDEX OF FIRST EQUATION TO BE SOLVED 
INDEX OF LAST EQUATION TO BE SOLVED 
FLAG FOR DUMP OPTIONS 
IL + 1 
IH - 1 

INPUT DEVICE 

FREE STREAM MACH NUMBER 

INDEX FOR MIXING LENGTH 

HAXIMUM NO. OF STORAGE LEVELS 

STORAGE LEVEL OF MIXING LENGTH IN 

GENERAL PURPOSE STORAGE 

INDEX FOR MACH NUMBER 


FORTRAN 

SYMBOL 

COMMON 

BLOCK 

DESCRIPTION 

MNL 

VAR 

STORAGE LEVEL OF MACH NO IN 



GENERAL PURPOSE STORAGE 

HOUT 

CIO 

OUTPUT DEVICE 

MREF 

REF 

REFERENCE MACH NUMBER 

MREEI 

REF 

INVERSE OF REFERENCE MACH NUMBER 

MSDD 

CIO 

TEMPORARY MASS STORAGE DEVICE 

MSD1 

CIOD 

GENERAL PURPOSE MASS STORAGE DEVICE 

MSD2 

CIOB 

GENERAL PURPOSE MASS STORAGE DEVICE 

MSGVARC25) 

PRNI 

TITLE OF VARIABLES TO BE PRINTED 

MU 

VAR 

INDEX FOR LAMINAR VISCOSITY 

MUINF 

FREE 

FREE STREAM LAMINAR VISCOSITY 

MUL 

VAR 

STORAGE LEVEL OF LAMINAR VISCOSITY IN 



GENERAL PURPOSE STORAGE 

MUREF 

REF 

REFERENCE VISCOSITY 

MUREFI 

REF 

INVERSE OF REFERENCE VISCOSITY 

MUT 

VAR 

INDEX FOR TURBULENT VISCOSITY 

MUTL 

VAR 

STORAGE LEVEL OF TURBULENT VISCOSITY IN 



GENERAL PURPOSE STORAGE 

MUINF 

FREE 

FREE STREAM MOLECULAR UEIGHT 

MUREF 

REF 

REFERENCE MOLECULAR UEIGHT 

MUREFI 

REF 

INVERSE OF REFERENCE MOLECULAR UEIGHT 

MZVAR 

PARAM 

MAXIMUM NUMBER OF SIORAGE VARIABLES 

NABC 

CCOM 

NABC = ID 

NANG 

METRIC 

ANGLE OF COORDINATE LINES RELATIVE TO HORIZONTAL 

NCPLD 

PARAM 

MAXIMIUN NUMBER OF COUPLED EQUATIONS 

NCPLD2 

CCOM 

NCPLDAA2 
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FORTRAN 

SYMBOL 

NCTR 

NCUP 

NDIFM 

NDIFMT 
NDIFM1 
NDIFP1 
NE( 2) 

NEQN(NEQS,2) 

NEQS 

NEY 

NEYM1 

NEZ 

NEZM1 

NE2S 

NFILE 

NGEOMV 

NH1 

NH12 

NH2 

NH21 

NH3 

NH3 1 

NH32 

NIIT 


COMMON DESCRIPTION 

BLOCK 

DIFCOM CENTER OF DIFFERENCE MOLECULE 

EON NUMBER OF COUPLED EQUATIONS 

PARAM MAX INIUM NUMBER OF GRID POINIS IN A DIFFERENCE 
MOLECULE 

DIFCOM 2ANDIFM 

DIFCOM NDIFM - 1 

DIFCOM NDIFM + 1 

GEOM NUMBER OF GRID POINTS IN Y AND Z DIRECTIONS 

EQN NUMBER OF COUPLED EQUATIONS TO BE SOLVED IN EACH 

ADI SWEEP 

PARAM MAXIMIUM NUMBER OF EQUATIONS IN CODE 
GRID NUMBER OF GRID POINTS IN Y DIRECTION 

GRID NEY - 1 

GRID NUMBER OE GRID POINTS IN Z DIRECTION 

GRID NEZ - 1 

CIO SENTINEL FOR SPREADING OF 2-D PROFILE TO 3-D 

REST SEQUENCE NUMBER OF RESTART INFORMATION 

METRIC NUMBER OE METRIC COEFFICIENTS AND DERIVATIVES 
METRIC INDEX FOR METRIC COEFFICIENT IN X-DIRECTION 
METRIC INDEX FOR DERIVATIVE OF X METRIC IN Y DIRECTION 
MEIRIC INDEX FOR METRIC COEFFICIENT IN Y-DIRECTION 
METRIC INDEX FOR DERIVATIVE OF Y METRIC IN X DIRECTION 
METRIC INDEX FOR METRIC COEFFICIENT IN Z-DIRECTION 
METRIC INDEX FOR DERIVATIVE OF Z METRIC IN X DIRECTION 

METRIC INDEX FOR DERIVATIVE OF Z METRIC IN Y DIRECTION 

CIO NO. OF FALSE MARCHING STEPS USED TO GENERATE 
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FORTRAN 

SYMBOL 

COMMON 

BLOCK 

DESCRIPTION 
THE INITIAL PROFILE 

NJD 

DIFCOM 

GRID POINT LOCATION FOR START OF SECOND SUEEP 

NMAXUF 

LAWW 

MAXIMIUM NUMBER OF ITERATIONS ALLOWABLE IN 
CALCULATIOON OF WALL SHEAR VELOCITY 

NN 

PARAM 

MAXIMIUM NUMBER OF GRID POINTS IN Y OR Z DIRECTION 

NPADI 

EQN 

NUMBER OF COUPLED AND UNCOUPLED EQUATIONS TO BE 
SOLVED DURING EACH ADI SWEEP 

NPTSX(4) 

INTBC 

NUMBER OF STREAMWISE LOCATIONS WHERE TRANSPIRATION 
DATA IS INPUT 

NPTSYZ<15,4) 

INTBC 

NUMBER OF CROSS-PLANE LOCATIONS WHERE TRANSPIRATION 
IS INPUT 

NPUNCH 

CIO 

PUNCH DEVICE 

NRGH2) 

DIFCQH 

NCTR POINTS FROM RIGHT OR TOP BOUNDARY 

NS 

GEOM 

LAST MARCHING STATION 

NSAVED 

REST 

SEQUENCE NUMBER OF RESTART STATIONS SAVED 

NSMAX 

PARAM 

2 GREATER THAN NS 

NUNERR 

CIO 

DEVICE FOR MGAUSS ERROR CHECK 

NVSOLV 

EQN 

NUMBER OF DEPENDENT VARIABLES 

NU0RD2(5O) 

STRAGE 

SIZE OF EACH COMMON BLOCK 

OMBWD 

LIN 

1.0 - BWD 

OMEGUF 

LAUU 

UNDER-RELAXAIION FACTOR FOR WALL FUNCTION 
FORMULATION 

P 

VAR 

INDEX FOR STATIC PRESSURE 

PCON1 

REF 

(GAMMA-1. 0)/GAMMA 

PC0N2 

REF 

0.5 A PC0N1 

PINF 

FREE 

FREE STREAM STATIC PRESSURE 
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FORTRAN 

SYMBOL 

COMMON 

BLOCK 

DESCRIPTION 

PL 

VAR 

STORAGE LEVEL OF STATIC PRESSURE IN 
GENERAL PURPOSE STORAGE 

PLTFLD<NN,NN f 8) 

ZPLOT 

GENERAL PURPOSE STORAGE FOR PLOT INFORMATION 

PR 

REF 

PRANDIL NUMBER 

PREE 

REF 

REFERENCE PRESSURE 

PREFI 

REF 

INVERSE OF REFERENCE PRESSURE 

PREPS 

vise 

PRANDIL NO. IN TURBULENT ENERGY DISSIPATION 
EQUATION 

PRESS ( 4 ) 

BOUND 

SPECIFIED PRESSURE AT BOUNDARY 

PRT 

REF 

TURBULENT PRANDIL NO. 

PRTKE 

vise 

PRANDIL NO. IN IURBOLENT KINECTIC ENERGY EQUATION 

PZERO 

FREE 

STAGNATION PRESSURE 

Q1D(8,NN) 

METRIC 

INTERMEDIATE STORAGE ARRAY FOR METRIC 
INFORMATION 

Q2D < 8 , NN ) 

METRIC 

INTERMEDIATE STORAGE ARRAY FOR METRIC 
INFORMATION 

R 

VAR 

INDEX FOR DENSITY 

RATLD 

LAUU 

EMPIRICAL NUMERICAL CONSTANT IN MIXING LENGIH 
COMPUTATION 

RE 

REF 

REYNOLDS NUMBER 

RE I 

REF 

INVERSE OF REYNOLDS NUMBER 

RE 12 

REF 

2.0 k RE I 

REPL 

FREE 

REYNOLDS NUMBER PER UNIT LENGTH 

RGAS 

REF 

GAS CONSTANT 

RHO(NDIFh) 

ADDR 

STORAGE OF DENSITY FOR FIRST DERIVATIVES 

RHOINF 

FREE 

FREE STREAM DENSITY 
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FORTRAN 

SYMBOL 

COMMON 

BLOCK 

DESCRIPTION 

RHOREE 

REF 

REFERENCE DENSITY 

RHOUL(NN,4> 

BOUND 

STORAGE FOR DENSITY ON THE BOUNDARY 

RHREEI 

REE 

INVERSE OF REFERENCE DENSITY 

RUNIV 

UNITS 

UNIVERSAL GAS CONSTANT 

SAVE(NEQS,NN) 

EON 

STORAGE FOR CHANGES DURING FIRST ADI SUEEP 

SN(NN) 

LIN 

STORAGE FOR N TH LEVEL TERMS 

SQ1 

ADD 

ADD CODE STREAMUISE LOCATION 

SQ2 

ADD 

ADD CODE STREAMUISE LOCATION 

SYSTEM 

ZPLOT 

SENTINEL FOR COORDINATE SYSTEM - PLOTS ONLY 

T 

VAR 

INDEX FOR STATIC TEMPERATURE 

TEMPS(NN f 4) 

BOUND 

STORAGE ARRAY FOR TEMPERATURE ON BOUNDARY 

TEMPSNINN, 4) 

BOUND 

STORAGE ARRAY FOR TEMPERATURE AT NTH STREAMUISE 



STATION - BOUNDARIES ONLY 

TINE 

FREE 

FREE STREAH STATIC TEMPERATURE 

TIILE(6) 

ZPLOT 

TITLE FOR PLOT FILE 

TKE 

VAR 

INDEX FOR TURBULENT KINETIC ENERGY 

IKEINE 

LAWU 

FREE STREAH TURBULENT KINETIC ENERGY 

TL 

VAR 

STORAGE LEVEL OF STATIC TEMPERATURE IN 



GENERAL PURPOSE STORAGE 

IREF 

REF 

REFERENCE IEMPERATURE 

TREEI 

REF 

INVERSE OF REFERENCE TEMPERATURE 

TT1 (2) 

GEOM 

MESH DISTRIBUTION FACTOR 

TT2(2) 

GEOM 

MESH DISTRIBUTION FACTOR 

TUALL(4) 

BOUND 

SPECIFIED IEMPERATURE AT BOUNDARY 

TUOD 

GEOM 

FLAG FOR TUO DIMENSIONAL LOGIC 

TZERO 

FREE 

STAGNATION TEMPERATURE 
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FORTRAN 

SYMBOL 


COMMON 

BLOCK 


DESCRIPTION 


12(2) 

GEOM 

U 

VAR 

UDUE ( 4 ) 

LAUU 

UDUM(NDIFM) 

ADDR 

UINF 

FREE 

UREF 

REF 

UREFI 

REF 

USCALE 

ADD 

USIAR<NN,4) 

LAUU 

UTIL(NDIFh) 

ADDR 

V 

VAR 

VELSQ(NN f NN) 


VKB 

LAUU 

VKC 

LAUU 

VNO( 15, 15,4) 

INIBC 

U 

VAR 

X(502) 

GEOM 

XENTR 

GEOM 

XG1 (NN,2) 

OPER 

XG2(NN, 2) 

OPER 

XOB ( 4 ) 

BOUND 

XVNO( 15,15,4) 

INTBC 


MESH DISTRIBUTION FACTOR 

INDEX FOR VELOCITY IN X-DIRECIION 

BOUNDARY LAYER THICKNESS SAMPLING CRITERIA 

STORAGE OF VELOCITY EOR FIRST DERIVATIVES 

FREE STREAM VELOCITY 

REFERENCE VELOCITY 

INVERSE OF REFERENCE VELOCITY 

METRIC SCALE FACTOR 

FRICTION VELOCITY ON SOLID WALL BOUNDARY 

STORAGE OF DEPENDENT VARIABLE FOR FIRST DERIVATIVES 

INDEX FOR VELOCITY IN Y-BIRECTION 

STORAGE FOR UAA2 + VAA2 + UAA2 - EQUIVALENCED TO 

PLTFLDU, 1,4) 

SECOND CONSTANT IN LOGARITHMIC LAU OF THE UALL 

VON KARMAN CONSTANT 

INPUT TRANSPIRATION RATES 

INDEX FOR VELOCITY IN Z-DIRECTION 

STREAMUISE LOCATION 

STARTING STREAMUISE LOCATION 

FIRST DERIVATIVES OF COMPUTATIONAL COORDINATES 
WITH RESPECT TO PHYSICAL COORDINATES 
SECOND DERIVATIVES OF COMPUTATIONAL 
COORDINATES WIIH RESPECT TO PHYSICAL COORDINATES 
INITIAL LOCATION OF SOLID OBSTACLE IN X DIRECTION 
STREAMUISE LOCATIONS UHERE TRANSPIRATION DATA IS 
INPUT 
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FORTRAN 

COMMON 

DESCRIPTION 

SYMBOL 

BLOCK 


XO (2) 

GEOM 

MESH DISTRIBUTION FACTOR 

Y(NN, NN,2) 


PHYSICAL DISTANCES FROM BOUNDARIES - EQUIVALENCED 
TO PLTFLDd, 1,2) 

YAW 

REF 

YAW ANGLE 

YPLUSL(NN,NN) 


STORAGE FOR RHO A UTAU / VISLAH - EQUIVALENCED TO 
PLTFLD( 1,1,1) 

YS(2,2) 

GEOM 

NONDIMENS IONAL EXTENTS OF COMPUTATIONAL DOMAIN 

YSAVE(NN, 2) 

GEOM 

COMPUTATIONAL COORDINATES 

YZPROF(NN) 

PRFILE 

TEMPORARY STORAGE ARRAY FOR PHYSICAL COORDINATES 

YZVN0U5, 15,4) 

INTBC 

CROSS-PLANE LOCATIONS WHERE TRANSPIRATION DATA IS 
INPUT 

ZNTRN(NDIN,3) 

ZEX1 

ARRAY FOR VIRTUAL MEMORY STORAGE 
NDIM = NN A NN A MZVAR A MLEVEL 

ZZ(M f L, K) 

VARZZ 

GENERAL PURPOSE STORAGE FOR DEPENDENT AND DERIVED 
VARIABLES - M = MZVAR, L = MLEVEL, K = NDIFM A NN 
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Sample Input and Sample Output 


Sample input and output for two- and three-dimensional cases are 

presented in Tables II - V, respectively. The two-dimensional input is that 

used in running the previously discussed Rose case (Ref. 36) while the 

three-dimensional input is that used in the running of the Bogdonoff case 

(Ref. 34). The Rose case input data is for an initial run (IRSTIN=0) with a 

restart to be written every 200 marching steps (IRSTOT). The case is to be 

run on our IBM virtual memory marching (IC0MP=3). The streamwise and 

transverse momentum as well as continuity and stagnation enthalpy (IHSTAG) 

version of the energy equation are to be solved. On the first boundary, the 

boundary type as well as the boundary condition are to change at streamwise 

distance 0.5618 (ASW( 1 )=0. 5618) . Wall temperatures on the 1 and 2 surfaces 

are 3.651111439 times the reference temperature (TWALL(1)=2*3. 651111439) . 

The reference length is 0.8666667 ft. (IUNITS*!) the Reynold's number per 

ft. is 5.22 x 10 6 , the free stream March number 3.88 and the free stream 

pressure is 5.9.8170527 lbt/ft 2 . Grid point packing about the 1 surface is 

to be moderate TTl=-0.80 while the packing around the 2 surface is to be 

considerably tighter TT2=0.95 99 grid points are utilized in the transverse 

direction (NE(1)=99) and an axisymmetric coordinate system generated by the 

ADD code (IGE0M=11) is to be utilized. The initial run is to be marched 200 

steps (NS=200) starting at a streamwise location of 0.2 (XENTR=0.2). An 

initial profile is to be supplied off the 2 surface. The boundary layer 

- 3 

thickness is 0.130 and the skin friction coefficient is 1.72 x 10 . A 

turbulent mixing length model is used (IMIXL— 1). Printout is given every 10 
steps ( IPRINT=10) and plot information is written every 2 steps (IPL0T=2). 

The output for the two-dimensional Rose case consists of NAMELIST information 
and geometric information and maximum change information and flowfield 
information at each tenth streamwise station. The NAMELIST information is 
provided as a means for the user to check the input data. The geometric 
information consists of the nond imensional computational distance (YSAVE) , 
the nondimensional physical distance Y from the lower surface, the three 
metrics and their derivatives and the X and Y physical location based on an 
(ADD code) absolute frame of reference. The maximum change information 
consists of the variable number (IVAR) the grid point position of the maximum 
change (JMAX and KMAX) and the value of the maximum change of the variable 
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from one marching station to the next. The next two pages consist of the 
flowfield information at a given strearawise location. All variables (except 
the pressure) are in a nondimens ional form with respect to the reference 
conditions which are displayed in NAMELIST LIST2 . The pressure terms are 
nondimensional ized with respect to the free stream pressure. The first page 
of flowfield output consists of the transverse grid point number, the 
computational position, the streamwise velocity component, the transverse 
direction velocity component, the density, the static or stagnation ( IHSTAG=0 
or IHSTAG=1 or 2) enthalpy, the turbulence kinetic energy, the static 
pressure, the static temperature, the effective viscosity and the Mach 
number. The second page of flowfield output again consists of the grid point 
number at the computational position followed by the subsonic (INDC=1) - 
supersonic (INDC=3) indicator, the stagnation temperature and pressure, the 
pressure coefficient, the laminar viscosity, the mixing length, the 
dissipation function and finally the cell Reynolds number. Following the 
flowfield information is the subsonic-supersonic grid point position 
indicator (ISS) with respect to the lower and upper surface and the boundary 
indicator (IBOUND). The ISS values tell the grid point where the flow 
transitions from subsonic to supersonic flow while the variable IBOUND tells 
the type of surface (JB0UND=1 corresponding to a wall and JB0UND=2 
corresponding to a nonwall). Finally, plot file information is displayed. 

The sample input for a three-dimensional restart case is presented in 
Table IV. Since this is a restart case the initial input values are retained 
as defaults and only variables that are going to be changed need appear. For 
this particular case, the results at step 350 (IRSTIN=350) are going to be 
marched 50 more steps to station 400 (NS=400). From this input stream it is 
possible to see how the streamwise step size was varied. The initial step 
size (which would have to be obtained from the initial input run stream) was 
decreased by 20% per step over the first 20 stations and allowed to reach a 
minimum step size of 0.002. At station 21 the step size was increased by 5% 
per step. From station 41 to 101 the step size remained constant. At 
station 101 the step size was again reduced by 20% per step with a minimum 
value of 0.0002. After station 251, the step size increased by 15% per step 
until it achieved a maximum step size of 0.005. A portion of the output for 
this case is presented in Table V. The initial three-dimensional output is 
the same as the two-dimensional output and is not presented here. The format 
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of the flowfield output is of a different form. Unlike the two-dimensional 
flowfield output, the three-dimensional output is controlled by the variable 
IVARPR. All of the variables that were printed inthe two-dimensional output 
can be obtained in three dimensions. For three dimensions, however, the 
output is in the form of a cross-sectional plane of output. The integer 
variables IY and IZ represent the transverse and spanwise grid point 
locations respectively while Z and Y are the corresponding compuational 
positions. Table IV is a portion of the output for the Bogdonoff case i.e., 
the cross plane distribution of the streamwise velocity (UVEL) and pressure 
(PRES). Other variables can (and were) printed out but for reasons of 
economy of space are not presented here. 
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FIGURE I- FORWARO MARCHING PROCEDURE - NOTATION 


9 



FIGURE 2- TYPICAL HIGH SPEED INLET 




primary ( x ) flow direction 
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FIGURE 3— TYPICAL COMPUTATIONAL DOMAIN FOR 
TWO DIMENSIONAL LAMINAR AND 
TURBULENT FLOW CALCULATION 
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Re/L = 10 /meter 
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Figure 4. Laminar Flat Plate - Streamwise Velocity Distribution. 




Laminar Flat Plate - Skin Friction Distribution 
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Figure 6, Supersonic Laminar Flat Plate - Streamwise Velocity Distribution, 
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Figure 7. Supersonic Laminiar Flat Plate - Skin Friction Distribution. 
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Figure 9. Supersonic Turbulent Flat Plate - Boundary Layer vs Re ( 
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(Reduced Navier-Stokes Analysis 



Figure 10. Supersonic Turbulent Flat Plate - Reduced Navier-Stokes C f vs Re 
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Figure 13* Hypersonic Corner Flow Schematic Representation of Corner Flow System. 
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Figure 14. Hypersonic Corner Flow - Computational Mesh in the Streamwise Plane. 
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Figure 17. Hypersonic Corner Flow Comparison of Skin Friction 
Coefficients on a Spanwise Surface. 
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Figure 18 . Hypersonic Corner Flow - Streamwise Development of the Stagnation 
Pressure Isobars. 
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Figure 22(a) - Gland 
Static 
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Figure 22(b) - Glancing Shock Wave - Boundary Layer Interaction 

Static Pressure at X = 7.6. 
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ave - Boundary Layer Interaction 




75 INCHES 



Figure 26(b) - Glancing Shock Wave - Boundary Layer Interaction 
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Figure 28 - Glancing Shock Wave - Boundary Layer Interaction 
Yaw Angle at X = 7.1. 
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Figure 30 - Nonorthogonal Channel Flow - Computational Mesh. 
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Figure 31 - Comparison of Computed Pressure Distribution 
on Upper Wall. 
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Figure 35 - Axisymmetric Inlet — • Mach 3.5 Inlet Model Schematic 



UJ 

u 

z 

2 

o 

< 

X 

< 


X 

w 

o 

X 

c 

o 


D 

D. 

g 

O 

U 


v£> 

m 

a) 

u 

3 

WJ 


1/A 


*30NVJLSia 1V10VB 



137 




isymmetric Inlet 








0.4 



141 


Figure AO - Ax i symmetric Inlet - Comparison of Cowl Static Pressure. 
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Figure 41 - Axisymmetric Inlet at Angle of Attack — Streanvwise Static 
Pressure Distribution on both Centerbody and Cowl on 
Windward Ray. 







Figure 43 - Axisymmetric Inlet Angle of Attack Streamwise Static Pressure 
Distribution on Both Centerbody and Cowl on Leeward Ray. 
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Figure 44 - Axisymmetric Inlet of Angle of Attack - Streamwise Development of 
Axial - Mach Number Profiles on Windward Ray. 
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Figure 45 - Axisymmetric Inlet of Angle of Attack - Streamwise Development of 
Mach Number Profiles on Waterline Ray. 




o 

tf> 


o 


o 


c 

<U 

a 

(X 

0 

t— I 

<u 

> 

<u 

Q 

d) 

CO 

•H 

1 

cd 

0) 

M 

4-» 

CO 



-I 

a • 

o 

\ 

cd 

to 

X 

4-t cd 
4J pS 


< 


•» 

-a 


tu 

CM V-l 

o cd 


o 

S 


z 

Q) Q) 


< 

H 

*— 4 <U 
60 hJ 

c 


V) 

< c 



o 


Q 

4-1 

cd co 

O 


0) 

cvi 

< 

4J »H 
<U 


X 

*— 1 vw 
C O 


<* 

M U 





vO 

a) 

3 

60 


o 

6 


147 




ocisymmetric Inlet of Ang] 
Actors in the Cross Sect: 








X/R L 

r/r l 

Slope 

Centerbody 

0.0 

0.0 

0.17633 

4.0 

0.70532 

0.17633 

4.1 

0.7228 


4.2 

0.7387 

0.144 

4.3 

0.7512 


4.4 

0.759 

0.052 

4.5 

0.7625 


4.55 

0.763 

0.0 

4.6 

0.7625 


4.65 

0.7611 


4.7 

0.7585 

-0.0646 

4.8 

0.7504 


4.9 

0.7391 

-0.1295 

5.1 

0.7120 


5.3 

0.6829 


5.5 

0.6525 

-0.153 

5.6 

0.6362 


5.7 

0.618 


5.8 

0.5973 


5.9 

0.5744 


6.0 

0.5467 


6.1 

0.5093 


6.2 

0.4564 


6.28 

0.4 

-0.794 

Cowl 

2.86 

1.0 

0.01745 

3. 1 

1.004188 

0.01745 

3.2 

1.0054 


3.4 

1.0051 

-0.011 

3.6 

0.99996 


3.8 

0.9882 


4.0 

0.9681 

-0.124 

4.1 

0.954 


4.2 

0.9364 

-0.1942 

4.25 

0.9261 


4.3 

0.9154 

-0.213 

4.4 

0.8949 


4.5 

0.8768 

-0.163 

4.55 

0.8695 


4.6 

0.864 

-0.093 

4.65 

0.86 


4.7 

0.8572 

-0.0485 

4.8 

0.8533 


4.9 

0.8511 













Table 1 - (Concluded) 


X/R L 

r/r l 

Slope 

Cowl 

5.0 

0.8502 


5.1 

0.85 

0.0 

5.6 

0.85 

0.0 

5.8 

0.8574 


5.9 

0.8646 


6.0 

0.8735 


6.1 

0.8839 

0.107 

6.2 

0.8946 


6.3 

0.9050 


6. A 

0.9145 


6.5 

0.9227 

0.0729 

6.6 

0.9299 


6.7 

0.9368 


6.8 

0.9435 


6.9 

0.95 

0.065 


Table 1 - Surface contours of the centerbody and cowl 
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TABLE II - Sample Two-Dimensional Input 


TOP RECORD 
****ROSE CASE**** 


1.0 
&REST 
ICOMP = 3, 

IRESTIN = 0, NFILE = 0, N SAVED = 0, 
IRSTOT = 200, 

JRSTOT - 17, 

SEND 
&LIST1 
IHSTAG = 1, 

IBOUND = 2, 1,2,2, 

IEQBC(l.l) = 2,2,13,13, 

IEQBC( 1 , 2) = 2,2,13,13, 

IEQBC( 1,3) = 2, 2, 2, 2, 

IEQBC( 1,4) = 16,16,16,16, 

IEQBC( 1,5) - 8, 8, 8, 8, 

JEQBC(l.l) = 44,21,11,11, 

JEQBC( 1,2) = 45,12,11,11, 

JEQBC( 1,3) = 12,12,2,2, 

JEQBC( 1,4) = 46,12,14,14, 

JEQBC( 1,5) = 47,17,17,17, 

ASW(l) = .5618, 

TWALL(l) = 2*3.651111439, 

&END 

&LIST2 

LREF = .08666667, 

RE PL = 5.22E+06, 

MINF = 3.88, 

PINF = 59.1870527 
PR = 0.710, 

IUNITS = 1, 

&END 

&LIST3 

XOB(l) = 0.5618, 

TWOD = .TRUE., 

TT1 = -0.80, 0.0, 

TT2 = 0.95, 0.0, 

YS( 1,1) - 1.0E-03, 

NE(1) = 99, 

IGEOM =11, 

DELX = 0.01, 

NS = 200, 

XENTR = 0.2, 

&END 

&LIST4 

IDIRP = 2.IBCP = 2, 

DELTAP(2) = 0.130, CFP(2) = 1.72E-03, 
IVISC = 3, 


TABLE II - Sample TWo-Dimensional Input (Continued) 


IPROF - 4, 

DELTAB(l) = 2*0.130, 
IMIXL = 1, 

&END 
&LIST5 
IPRINT = 10, 

IPLOT = 2, 

&END 


Table III - Sample Two-Dimensional Output 
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TABLE IV - Sample Three-Dimensional Input 


TOP RECORD 
BOG TEST CASE 


1.0 
&REST 
ICOMP = 3, 

IRSTIN = 350, 

IRSTOT = 50, 

JRSTOT = 17, 

NFILE = 1, 

NSAVED = 0, 

&END 

&LIST1 

SEND 

&LIST2 

&END 

&LIST3 

NS = 400, 

IAP = 1,21,41,101,251,5*10000, 

AP = 0.8,1.05,1.0,0.8,1.15,5*1.0, 

DXMIN = 0.002,0.0,0.0,0.002,6*0.0, 

DXMAX = 4*100000.0,0.005,5*100000.0, 

&END 

&LIST4 

&END 

&LIST5 

&END 
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